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1.0  Executive  Summary 

This  project  focused  on  improving  the  realism  of  combustion  chemistry  models  used  in 
existing  state-of-the-art  computational  fluid  dynamics  (CFD)  software  that  is  currently 
used  in  the  design  of  scramjet  combustors.  The  result  is  a  modeling  tool  that  is  a 
significant  improvement  over  existing  capabilities  and  the  ability  to  create  chemistry 
models  suitable  for  use  in  CFD  from  published  detailed  chemical  kinetic  mechanisms. 

Capabilities  were  developed  and  demonstrated  for  1)  creating  accurate  reduced  chemical 
kinetic  mechanisms  for  hydrocarbon  fuels  of  interest  for  high  speed  propulsion,  2) 
integrating  the  reduced  mechanisms  with  in-situ  adaptive  tabulation  (ISAT)  and  the  eddy 
dissipation  concept  (EDC)  turbulence-chemistry  interaction  model  into  the  high-speed 
CFD  code  VULCAN  (Viscous  Upwind  ALgorithm  for  Complex  Flow  ANalysis),  and  3) 
simulating  combusting  flows  with  complex  aviation  fuel  chemistry  in  scramjet  combustor 
geometries. 

The  key  results  from  the  Phase  II  work  included  the  following: 

Development  of  accurate  reduced  mechanisms  for  combustion  kinetics  of  hydrogen, 
ethylene,  « -heptane,  and  a  six-specie  surrogate  blend  representing  JP-8. 

Improvements  to  the  CARM  (Computer  Aided  Reduction  Mechanism)  software  for 
automated  reduction  of  chemical  kinetic  mechanisms,  resulting  in  improved 
robustness  and  decreased  computational  times. 

Implementation  of  four  reduced  mechanisms  into  VULCAN. 

Implementation  of  ISAT  chemical  rate  tabulation  technique  into  VULCAN. 
Implementation  of  EDC  turbulence  model  (with  ISAT)  into  VULCAN. 

Completion  of  scramjet  combustor  simulations  demonstrating  improved  chemical 
descriptions  from  reduced  mechanisms  compared  to  previously  used  global  chemistry 
models. 

ISAT  results  showing  speed  increases  up  to  a  factor  of  36  for  the  ethylene  case. 
Two-dimensional  combustion  simulations  in  VULCAN  demonstrating  mechanism 
reduction  and  implementation  for  hydrogen,  ethylene,  « -heptane,  and  JP-8  (up  to  35 
species  in  the  reduced  mechanism). 

Three-dimensional  scramjet  combustor  simulations  run  for  ethylene  and  « -heptane 
(up  to  20  species  in  the  reduced  mechanism),  demonstrating  that  reduced  mechanisms 
and  ISAT  can  be  used  for  CFD  simulations  of  realistic  3-  dimensional  scramjet 
geometries  with  hydrogen  fuels. 

Overall,  the  research  completed  during  this  project  has  demonstrated  newly  developed 
mechanism  reduction  techniques  applied  to  newly  acquired  knowledge  of  hydrocarbon 
combustion  kinetics.  These  techniques  were  integrated  into  a  practical  engineering  tool 
to  produce  a  method  for  calculation  of  realistic  combustion  chemistry  in  high-speed  CFD. 
This  capability  allows  the  U.S.  Air  Force  to  better  address  the  needs  of  its  high-  speed 
propulsion  programs,  specifically  the  simulation  of  hydrocarbon  fuels  in  a  scramjet 
combustor. 
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2.0  Introduction 

2.1  Purpose  of  Work 

Computational  modeling  has  become  a  central  analytical  tool  used  in  the  design  of 
practical  combustion  systems  for  increased  performance,  improved  efficiency,  and  low 
emissions.  CFD  simulation  is  particularly  important  in  the  aerospace  industry  in  reducing 
prototype  and  testing  costs  and  the  time  needed  to  bring  products  to  market.  Accurate 
simulation  of  chemical  processes  in  aircraft  engines  is  critical  because  chemistry 
determines  the  rate  of  heat  release,  flame  stability,  ignition,  extinction,  and  pollutant 
emissions.  However,  fully  detailed  chemical  kinetic  descriptions  of  hydrocarbon 
oxidation  may  require  the  tracking  of  hundreds  of  chemical  species  and  thousands  of 
reaction  steps.  CPU  and  memory  limitations  prohibit  implementation  of  full  detailed 
chemistry  of  practical  fuels  into  3-D  CFD  simulations,  even  using  the  latest  massively 
parallel  computers. 

Improved  methods  for  computing  combustion  in  high  speed  flows  are  of  particular 
interest  to  the  U.S.  Air  Force.  In  particular,  the  U.S.  Air  Force  HyTech  program  (Mercier 
and  Ronald,  1996)  is  one  example  of  the  need  for  a  reliable  representation  of  reduced 
kinetic  models  within  CFD  simulations.  The  goal  of  the  HyTech  program  is  to  develop 
and  demonstrate  a  hydrocarbon-fueled  engine  capable  of  accelerating  an  air-launched 
missile  to  Mach  8.  The  development  of  a  successful  scramjet  combustor  has  proven  to  be 
difficult  due  in  large  part  to  structural  limitations  caused  by  high  heat  transfer  rates,  and 
short  residence  times  for  fuel  and  air  mixing  (Heiser  and  Pratt,  1994,  and  Curran  et  al., 
1996).  CFD  is  routinely  used  to  analyze  the  combustor  geometry  concepts  for  these 
vehicles  both  before  and  after  tests  (Baurle  et  al.,  1997,  1998a,  1998b,  2000).  A  large 
portion  of  the  design  phase  of  these  combustors  must  rely  heavily  on  CFD  because  of 
difficulties  in  measuring  all  of  the  relevant  flow  parameters  that  affect  scramjet 
combustor  performance. 

The  grid  requirements  for  solving  the  Reynolds-averaged  Navier-Stokes  equations  on 
scramjet  combustor  geometries  can  easily  exceed  several  million  grid  points.  This 
limitation  places  a  significant  constraint  on  the  size  of  the  chemistry  models  that  can  be 
considered.  Calculations  have  shown  that  existing  global  mechanisms,  even  for  simple 
hydrocarbon  fuels  such  as  ethylene,  are  woefully  inadequate  (Baurle  et  al.,  1998b).  Thus, 
previous  CFD  modeling  of  scramjet  combustion  has  been  very  limited  in  terms  of 
including  adequate  chemical  representations  of  the  combusting  flow.  Further,  the  effect 
of  turbulence  on  chemical  reactions  has  not  been  adequately  addressed.  The  inclusion  of 
turbulence-chemistry  interactions  in  simulations  further  exacerbates  the  CPU  and 
memory  limitations. 

Recent  advances  in  chemical  and  numerical  techniques  have  helped  to  address  some  of 
these  simulation  challenges.  Detailed  chemical  kinetic  descriptions  of  the  oxidation  of 
hydrocarbon  blends  representing  practical  fuels  are  rapidly  improving  (Dagaut,  1994; 
Maurice,  1996;  Curran  et  al.,  1998;  and  Violi  et  al.,  2002).  These  detailed  mechanisms, 
though  impractical  for  implementation  into  a  CFD  simulation,  are  the  starting  point  for 
creating  reduced  chemical  kinetic  mechanisms  which  can  be  implemented  into  CFD 
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simulations.  Furthermore,  techniques  are  now  available  to  automatically  create  reduced 
mechanisms  that  approximate  the  results  of  detailed  chemical  kinetic  descriptions  using 
many  fewer  species,  and  thus  less  CPU  time  and  memory.  These  advances  have  the 
potential  to  be  combined  with  new  numerical  techniques,  such  as  ISAT,  to  further  reduce 
CPU  time  and  thus  promote  tractable  solutions  for  3-D,  high-speed,  chemically  reacting 
flows. 

This  program  was  designed  specifically  to  utilize  recent  advances  in  chemical  kinetic 
descriptions  and  numerical  techniques  to  enhance  modeling  capabilities  for  scramjet 
engines  firing  hydrocarbon  jet  fuels.  Beyond  this,  the  capabilities  developed  in  this 
program  can  also  be  applied  to  other  aspects  of  the  U.S.  Air  Force’s  high-speed 
propulsion  programs  as  well  as  to  CFD  modeling  of  combusting  flows  in  general.  Two 
examples  of  additional  applications  that  could  benefit  from  program  results  are  ignition 
delay  studies  for  supersonic  flows  and  prediction  of  atmospheric  pollutants  from  fossil 
fuel  combustion. 

A  more  general  area  where  the  need  for  improved  representations  of  finite-rate  chemistry 
is  particularly  evident  is  in  predicting  the  emission  of  atmospheric  pollutants,  such  as 
oxides  of  nitrogen  (NOx),  from  the  combustion  of  fossil  fuels.  Identifying  methods  for 
reducing  NOx  emissions  from  combustion  devices  such  as  gas  turbines,  internal 
combustion  engines,  fossil  fuel- fired  utility  boilers,  and  industrial  process  heaters  has 
become  tremendously  important  to  U.S.  industry  in  an  effort  to  meet  continually 
tightening  air  quality  standards.  Integration  of  more  accurate  representations  of  chemistry 
into  specialized  CFD  tools,  based  on  information  gained  from  this  research,  helps 
improve  the  ability  to  solve  complex  industrial  combustion  problems  and  related 
pollutant  emission  concerns. 

2.2  Related  Research 

This  SB1R  Phase  11  work  builds  on  research  previously  completed  in  Phase  1  as  well  as 
research  performed  in  parallel  programs  at  Reaction  Engineering  International  (RE1) 
funded  by  the  U.S.  Army,  National  Science  Foundation  (NSF),  and  U.S  .industry. 

The  historical  reasons  that  CFD  codes  have  typically  incorporated  only  grossly  simplified 
descriptions  of  the  chemistry,  as  outlined  in  the  previous  section,  are  not  as  valid  today  as 
they  were  in  the  past.  Advancements  have  been  made  in  computer  software,  numerical 
methods  and  hardware  that  provide  a  framework  for  increasingly  detailed  chemistry  to  be 
implemented  into  3-D  CFD  simulations  to  address  practical  engineering  problems. 
Validated,  detailed  mechanisms  for  the  large  hydrocarbons  that  comprise  aviation  fuels 
have  been  available  for  some  time  now  (Dagaut,  1994;  Maurice,  1996;  Curran  et  al., 

1998;  and  Violi  et  al.,  2002).  These  mechanisms  represent  the  state  of  the  art  in  terms  of  a 
chemical  kinetic  description  of  the  combustion  of  large  hydrocarbons.  However,  use  of 
fully  detailed  mechanisms  for  hydrocarbon  fuels  in  such  simulations  is  still  unreachable 
for  the  foreseeable  future.  Great  strides  have  been  made  that  improve  our  capability  to 
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rapidly  develop  reduced  descriptions  of  detailed  chemical  mechanisms  while  retaining 
the  necessary  fidelity  to  the  detailed  chemistry. 

The  Phase  I  SBIR  program  successfully  demonstrated  that  a  chemical  kinetic  mechanism 
for  a  large  hydrocarbon  fuel,  such  as  those  that  comprise  aviation  fuels,  can  be  quickly 
reduced  to  a  size  amenable  for  implementation  in  CFD  codes  using  automated  techniques 
such  as  CARM  (Chen,  1997).  The  success  of  the  reduction  strategy  in  the  Phase  I 
program  was  demonstrated  by 

Successful  application  to  a  range  of  large  detailed  mechanisms  for  practical  fuels 
Comparison  to  experiments 

Comparison  to  results  of  full  mechanisms  for  simple,  spatially  homogeneous 
combustion  configurations  such  as  perfectly  stirred  reactors  (PSRs)  and  plug  flow 
reactors  (PFRs) 

The  ability  of  the  reduced  mechanism  to  capture  heat  release,  stability,  and  ignition 
delay  behavior. 

During  the  Phase  I  program,  which  concluded  in  March  1999,  REI  successfully  produced 
reduced  mechanisms  with  various  numbers  of  species  and  reaction  steps  for  ethylene  and 
from  two  different  detailed  mechanisms  for  w-heptane,  a  large  hydrocarbon  fuel  similar  to 
aviation  fuel  constituents.  More  complex  fuels  (heavier  alkanes  and  aromatics)  require 
more  species  to  accurately  represent  their  chemistry.  However,  since  large  hydrocarbons 
break  down  rapidly  into  smaller  pieces  under  combustion  conditions,  the  number  of 
scalars  needed  for  accurate  description  of  the  reaction  kinetics  can  be  expected  to 
increase  slowly  with  increasing  fuel  molecular  weight. 

Beyond  the  Phase  I  program,  REI  has  led  research  projects  in  several  related  areas 
between  the  end  of  the  Phase  I  effort  and  the  present.  These  include  the  following: 

A  U.S.  Department  of  Defense  (DoD)  Army  Phase  II  SBIR  for  the  development  of 
detailed  and  reduced  mechanisms  and  application  in  CFD  codes  for  problems  of 
interest  to  the  DoD  (Denison  et  al.,  2003) 

U.S.  industry-funded  development  and  application  of  reduced  mechanisms  for 
prediction  of  NOx  reduction  in  fossil-fuel  fired  utility  boilers  (Cremer  et  al.,  2000a, 
2000b, 2001) 

NSF  Phase  I  and  II  SBIRs  for  investigation  of  advanced  Newton-Krylov  solver 
techniques  to  reduce  computational  times  in  chemically  reacting  flow  simulations, 
especially  those  using  reduced  mechanisms  created  using  CARM  (Bockelie  et  al., 
2003) 

A  U.S.  DoD  Air  Force  Phase  I  SBIR  for  numerical  and  experimental  investigation  of 
soot-reducing  fuel  additives  in  JP-8  jet  fuel  combustion,  including  implementation  of 
reduced  JP-8  kinetics  and  soot  models  into  the  UNICORN  CFD  code  (Montgomery  et 
al.,  2003a) 

NSF  Phase  I  and  II  SBIRs  for  development  of  a  computer-assisted  reduction 
mechanism  problem  solving  environment  (CARM-PSE)  that  allows  easy  and 
automated  comparison  of  detailed  and  reduced  chemistry  over  a  multidimensional 
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parameter  space  of  interest  (Montgomery  et  al.,  2002a).  During  Phase  II,  automated 
reduced  mechanism  optimization  capability  was  added  and  optimized  reduced 
mechanisms  integrated  into  the  commercial  CFD  code,  Fluent,  and  the  internal 
combustion  engine  research  code,  KIVA  . 

Information  and  techniques  developed  in  these  programs  has  helped  further  the  research 
on  this  Phase  II  program,  particularly  with  respect  to  the  development  and  application  of 
reduced  mechanisms  in  chemically  reacting  flows. 

2.3  Phase  II  Objectives 

The  overall  goal  for  this  SBIR  Phase  II  program  was  to  provide  the  U.S.  Air  Force  with 
an  improved  computational  tool  for  simulating  jet  fuel  combustion  in  a  scramjet 
combustor.  This  goal  was  met  by  achieving  the  following  objectives: 

Identify  detailed  mechanisms  best  describing  oxidation  of  hydrocarbon  fuels  for 
scramjets 

Produce  reduced  mechanisms  for  the  detailed  hydrocarbon  mechanisms  and  verify 
reduced  mechanism  accuracy  by  comparison  with  detailed  mechanism  results 
Implement  the  reduced  mechanisms  into  a  CFD  code  suitable  for  performing  2-D  and 
3-D  scramjet  combustor  simulations 

Validate  the  implementation  against  available  experimental  data  and  with  detailed 
chemistry  in  simplified  geometries 

Implement  and  test  new  chemistry  tabulation  techniques  to  reduce  computation  times 
in  the  CFD  code 

Evaluate  the  impact  of  turbulence-chemistry  interactions  on  predictions. 

2.4  Key  Accomplishments  of  Phase  II  Work 

The  Phase  II  objectives  were  addressed  by  focusing  on  the  following  five  areas  of 
research  and  development: 

Identification  of  detailed  chemical  kinetics  mechanisms  for  simple  and  complex 
hydrocarbon  fuels  appropriate  for  use  in  simulating  scramjet  combustion 
Formulation  of  reduced  chemical  kinetic  mechanisms  that  represent  the  detailed 
kinetic  behavior  of  hydrocarbon  fuels  using  many  fewer  species  than  the  detailed 
mechanism 

Implementation  of  a  chemical  rate  tabulation  technique  (ISAT)  that  significantly 
reduces  CPU  time  when  computing  chemical  reactions  with  reduced  mechanisms 
Implementation  of  reduced  mechanisms  and  tabulation  techniques  into  the  VULCAN 
CFD  code 

Preliminary  evaluation  of  turbulence-chemistry  impacts  on  VULCAN  simulations 
based  on  the  EDC  modeling  and  a  probability  density  function  (PDF)  post-processing 
technique. 

The  Phase  II  objectives  were  accomplished  by  integrating  key  technologies,  developed 
and/or  implemented  as  part  of  the  process  outlined  above,  to  develop  a  capability  for 
modeling  turbulent  reactive  flows  that  has  not,  as  yet,  been  fully  realized  for  high-speed 
reacting  flow  applications.  The  improved  simulation  capabilities  were  evaluated  using  2- 
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D  simulations  to  verify  implementation  and  test  result  sensitivities  and  3-D  simulations  to 
demonstrate  predictions  for  a  scramjet  combustor. 

Key  results  from  the  Phase  II  work  included  the  following: 

Development  of  reduced  mechanisms  for  combustion  kinetics  of  hydrogen,  ethylene, 
/7-heptane,  and  a  six-specie  surrogate  blend  representing  JP-8 
Comparison  with  detailed  results  showed  excellent  agreement  over  a  wide  range  of 
temperatures  and  equivalence  ratios 

Improvements  to  CARM  software  for  automated  reduction  of  chemical  kinetic 
mechanisms  resulting  in  improved  robustness  and  decreased  computational  times 
Implementation  of  four  reduced  mechanisms  into  VULCAN 
Implementation  of  ISAT  chemical  rate  tabulation  technique  into  VULCAN 
Implementation  of  EDC  turbulence  model  (with  ISAT)  into  VULCAN 
Utilization  of  PDF  post-processing  technique  for  estimation  of  turbulence-chemistry 
interaction 

Simulation  results  showing  chemical  kinetic  effects  missing  from  simpler  global 
chemistry  models 

ISAT  timing  results  showing  speed  increases  up  to  a  factor  of  36 
Turbulence-chemistry  interactions  showing  significant  impact  on  the  flame  front; 

2- D  combustion  simulations  run  for  hydrogen,  ethylene,  n-heptane  and  JP-8  (up  to  35 
species  in  the  reduced  mechanism) 

3- D  scramjet  combustor  simulations  run  for  ethylene  and  n-heptane  (up  to  20  species 
in  the  reduced  mechanism). 

This  research  has  demonstrated  the  application  of  newly  developed  mechanism  reduction 
and  rate  tabulation  techniques  to  newly  developed  hydrocarbon  combustion  kinetics  and 
the  subsequent  integration  of  these  techniques  into  a  practical  engineering  tool  to  produce 
a  method  for  calculation  of  realistic  chemistry  in  CFD.  This  capability  allows  the  U.S. 

Air  Force  to  better  address  the  needs  of  its  high-speed  propulsion  programs,  specifically 
the  simulation  of  hydrocarbon  fuels  in  a  scramjet  combustor. 

3.0  Technical  Approach 

This  section  provides  an  overview  of  the  technical  approach  and  key  technical  tasks  in 
this  research.  Figure  1  illustrates  conceptually  the  technical  approach  and  task 
dependencies  for  this  program.  The  technical  research  and  code  development  were 
divided  into  five  development  tasks  plus  the  evaluation  of  2-D  and  3-D  combustion 
simulation  results. 
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Figure  1.  Technical  Approach  Used  in  Program. 

The  five  development  tasks  were  as  follows: 

1 .  Selection  of  detailed  chemical  kinetic  mechanisms  for  fuels  of  interest 

2.  Formulation  of  reduced  chemical  kinetics  mechanisms  using  CARM 

3.  Implementation  of  IS  AT  routines  for  use  with  reduced  mechanisms  and  VULCAN 

4.  Implementation  of  EDC  routines  for  use  with  reduced  mechanisms  and  VULCAN 

5.  Implementation  of  reduced  mechanisms,  IS  AT  and  EDC  into  VULCAN. 

Once  the  development  and  implementation  tasks  were  completed,  VULCAN 
simulations  were  run  to  evaluate  the  performance  of  the  reduced  mechanism, 

ISAT,  and  EDC  routines,  respectively.  The  following  sections  provide  a  summary  of 
each  of  the  development,  implementation,  and  simulation  tasks  for  this  program. 

3.1  Detailed  Mechanism  Selection 

This  task  involved  the  selection  of  detailed  chemical  kinetic  mechanisms  for  simple  and 
complex  aviation  fuels.  Reduced  mechanisms  can  be  no  more  accurate  than  the  detailed 
chemical  kinetic  models  on  which  they  are  based.  Thus,  accurate  description  of  detailed 
chemical  kinetic  mechanisms  is  a  necessary  first  step  for  each  fuel  of  interest. 
Identification  of  detailed  mechanisms  was  divided  into  two  categories,  simple  fuels 
(hydrogen  and  ethylene)  and  complex  fuels  (heptane  and  JP-8).  Section  4  of  this  report 
describes  the  selection  of  detailed  mechanisms  for  reduction  and  eventual  incorporation 
into  VULCAN.  Detailed  mechanisms  were  taken  from  existing  literature.  Some 
modifications  for  compatibility  with  CARM  software  were  required  depending  on 
original  detailed  mechanism  format. 

3.2  Formulation  of  Reduced  Mechanisms 

The  detailed  chemical  kinetic  mechanisms  identified  in  the  previous  task  were  utilized  to 
formulate  reduced  chemical  kinetic  mechanisms  using  CARM.  This  allows  a  great  deal 
of  chemical  kinetic  information  contained  in  the  detailed  mechanisms  to  be  replicated  in 
a  computationally  tractable  form.  Reduced  mechanisms  were  developed  for  simple 
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(hydrogen,  ethylene)  and  complex  (heptane  and  JP-8)  hydrocarbon  fuels.  Accuracy  of  the 
reduced  mechanisms  was  validated  by  comparison  to  detailed  chemistry  over  a  range  of 
conditions  and,  where  available,  experimental  data.  Details  of  the  formulation  approach 
and  results  are  given  in  Section  5  of  this  report. 

3.3  Implementation  of  ISAT  Routines 

The  ISAT  technique  provides  a  method  for  accessing  chemical  properties  without 
resorting  to  direct  calculation  of  reaction  rates  or  direct  tabulation  of  properties.  This 
results  in  significant  CPU  savings.  The  ISAT  subroutines,  developed  originally  for  a  PDF 
method  (Chen  et  al.,  2000),  were  modified  and  used  to  tabulate  the  chemical  reaction 
rates  from  both  the  CARM-produced  reduced  mechanisms  and  the  EDC  model. 
Background  on  the  ISAT  approach  is  given  in  Section  6  of  this  report. 

Since  this  was  the  first  use  of  ISAT  for  this  type  of  CFD  application,  numerous  timing 
and  accuracy  studies  were  conducted  for  different  ISAT  settings.  VULCAN  timing 
studies  were  performed  with  and  without  ISAT  for  many  of  the  reduced  mechanisms  in 
order  to  develop  a  sense  of  appropriate  accuracy  and  tolerance  settings  to  be  used  within 
ISAT. 


3.4  Turbulence-Chemistry  Interactions 

The  modeling  of  the  mean  chemical  source  term  has  often  been  considered  to  be  the  main 
problem  in  CFD  simulations  of  turbulent  combustion  (Peters,  2000).  Turbulent 
combustion  models  can  be  placed  into  one  of  following  categories:  eddy-breakup  (EBU) 
and  eddy-dissipation  models,  PDF  model  (including  assumed  PDF  method  and  transport 
PDF  method),  laminar  flamelet  models,  conditional  moment  closure,  and  the  linear  eddy 
model  (LEM).  Transport  PDF  and  LEM  models  are  accurate,  but  they  are  also  extremely 
expensive.  Application  of  the  laminar  flamelet  model  is  limited  when  more  than  one  fuel 
stream  is  present.  The  CPU  cost  of  the  assumed  PDF  model  increases  dramatically  when 
applied  to  the  chemical  mechanisms  without  an  analytical  reaction  rate  expression,  such 
as  the  CARM-created  reduced  mechanisms.  In  this  project,  the  turbulence-chemistry 
interactions  are  assessed  via  two  approaches. 

First,  an  EDC  model  was  implemented  into  VULCAN  and  configured  to  work  with 
ISAT.  This  was  used  to  test  the  sensitivity  of  VULCAN  predictions  to  turbulence- 
chemistry  interactions  within  the  limitations  of  the  EDC  model  (see  Section  7.1).  Second, 
a  PDF  post-processor  technique  was  applied  to  VULCAN  results.  This  technique  is 
described  in  Section  7.2.  The  PDF  post-processor  code  reads  the  flow  field  from  the 
VULCAN  results  and  recalculates  the  turbulent  species  concentrations  and  temperature 
field  with  a  PDF  model,  while  the  velocity  and  density  are  not  updated.  The  results 
obtained  with  the  EDC  and  PDF  models  are  then  compared  with  those  from  laminar 
chemistry  to  estimate  the  impact  of  turbulence  effects  on  flow  properties  and  CPU  times. 

3.5  Implementation  into  VULCAN 

The  simulation  software  utilized  in  this  research  for  modeling  scramjet  combustion  was 
VULCAN.  VULCAN  is  a  time-dependent,  compressible,  Navier-Stokes  flow  solver  that 
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is  currently  maintained  and  distributed  (restricted  to  the  U.S.)  by  the  Hypersonic 
Airbreathing  Propulsion  Branch  of  NASA  Langley  Research  Center.  VULCAN  is  a 
multigrid,  multiblock,  structured,  finite-volume  code,  developed  for  solving  the  spatially 
elliptic  and  parabolized  forms  of  the  equations  governing  3-D,  turbulent,  calorically 
perfect  and  nonequilibrium  chemically  reacting  flows.  The  reduced  chemical 
mechanisms,  ISAT  chemical  property  routines,  and  EDC  turbulence  model  routines 
developed  for  this  research  were  implemented  into  VULCAN  (see  VULCAN  web  page). 

Implementation  of  reduced  mechanisms  in  the  VULCAN  code  focused  in  two  areas. 

First,  a  numerical  Jacobian  capability  was  implemented  in  VULCAN.  This  was  necessary 
to  allow  VULCAN  to  use  the  (nonanalytical)  chemical  rate  source  terms  from  CARM- 
formulated  reduced  mechanisms.  Second,  the  reduced  mechanisms  for  various 
hydrocarbon  fuels  were  implemented  and  tested  for  numerical  accuracy  and  robustness. 

The  ISAT  routines  functioned  as  an  interconnection  between  VULCAN  and  the  reduced 
mechanism  routine  from  CARM.  The  ISAT  subroutine  provided  not  only  the  reaction 
rate  vector,  but  also  the  Jacobian  matrix  required  by  VULCAN.  This  provided  an 
opportunity  for  significant  CPU  savings  over  direct  evaluation  of  the  reaction  rates. 

Implementation  of  the  EDC  model  into  VULCAN  was  very  similar  to  the  implementation 
of  the  reduced  mechanism  with  ISAT.  In  this  formulation,  the  target  functions  of  ISAT 
were  the  turbulent  mean  reaction  rates,  and  the  independent  variables  for  ISAT  were  the 
species  mass  fractions,  pressure,  temperature  and  the  life  time  of  the  fine-scale  structures. 
This  required  only  one  more  variable  than  that  of  ISAT  with  laminar  chemistry.  This 
allowed  manageable  CPU  times  for  the  EDC  calculations. 

The  PDF  post-processor  utilized  information  from  the  VULCAN  flow  field  results  to 
calculate  turbulence-chemistry  effects  as  a  post-processor.  Thus  there  was  no  VULCAN 
implementation  required. 

3.6  2-D  and  3-D  Simulations 

2-D  VULCAN  simulations  were  used  to  verify  correct  implementation  of  the  reduced 
mechanism,  ISAT  and  EDC  routines  in  VULCAN,  and  to  test  the  sensitivity  of  predicted 
results  (accuracy,  CPU  time)  to  variations  in  new  routine  parameters  such  as  search 
tolerances  and  numerical  convergence.  Most  of  the  project  simulations  completed  were  2- 
D.  Once  the  correct  implementation  and  usage  of  the  new  routines  in  VULCAN  were 
established  with  the  2-D  simulations,  a  few  3-D  simulations  were  performed  to 
demonstrate  the  capabilities  of  the  newly  combined  technologies.  3-D  cases  were  run  for 
ethylene  and  n -heptane  fuels  in  a  scramjet  cavity  flameholder. 

4.0  Selection  of  Detailed  Mechanisms 

Reduced  mechanisms  can  be  no  better  than  the  detailed  chemical  kinetic  models  on 
which  they  are  based.  This  section  describes  our  selection  of  detailed  mechanisms  for 
reduction  and  eventual  incorporation  into  VULCAN. 
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4.1  Detailed  Mechanisms  for  Simple  Fuels 

4.1.1  Hydrogen 

The  reduced  mechanism  was  based  on  H2-O2  reactions  extracted  from  the  large 
hydrocarbon  mechanism  of  Curran  et  al.  (1998).  The  hydrogen-oxygen  (H-O)  reactions 
in  most  modem  detailed  mechanisms  differ  little  and  almost  any  detailed  mechanism 
could  have  been  chosen.  The  Curran  mechanism  is  a  widely  used  model  of  large 
hydrocarbon  combustion,  and  the  H-0  subsystem  is  believed  to  be  quite  reliable.  The 
H-0  reactions  are  modeled  using  nine  species  and  27  reactions. 

4.1.2  Ethylene 

Calculations  were  performed  to  compare  the  ability  of  two  detailed  chemical  kinetic 
mechanisms  (Marinov  et  al.,  1998,  and  Wang  et  al.,  1999)  to  compute  ignition  delay  of 
ethylene-oxygen  mixtures  as  compared  to  experiments  (Baker  and  Skinner,  1972,  and 
Colket  and  Spadaccini,  1999).  Three  of  the  cases  measured  by  Baker  and  Skinner  were 
modeled.  The  initial  pressure  in  all  of  these  cases  was  3  atm.  The  gas  mixture  for  these 
cases  consisted  of  1  percent  (by  moles)  ethylene  with  3,  1.5,  or  6  percent  O2,  and  the 
balance  argon.  These  cases  correspond  to  E.R.s  of  1.0,  2.0,  and  0.5,  respectively.  Baker 
and  Skinner  report  only  Arrhenius-form  expressions  fit  to  their  data  so  individual  data 
points  can’t  be  plotted.  Figure  2  compares  Baker  and  Skinner’s  measurements  with 
calculations  using  the  two  detailed  mechanisms  and  the  20-specie  reduced  mechanism 
described  in  Section  5.2.2.  Three  cases  measured  by  Colket  and  Spadaccini  were  also 
compared  to  predictions  of  the  detailed  chemical  kinetic  mechanisms  (not  shown).  These 
measurements  have  0.7  percent  ethylene  with  O2  mole  fractions  of  4.2,  2.8,  and  2.1 
percent  with  the  balance  argon.  This  gives  E.R.s  of  0.5,  0.75,  and  1.0,  respectively.  The 
initial  pressures  range  from  5.27  to  7.89  atm.  Similar  levels  of  agreement  for  the  detailed 
and  reduced  mechanisms  were  found  for  the  Baker  and  Skinner  measurements. 

Both  detailed  mechanisms  give  reasonable  agreement  to  the  measurements.  However, 
the  Wang  mechanism  is  noticeably  superior  for  all  of  the  measurements  simulated.  The 
Wang  mechanism  also  has  the  advantage  of  containing  fewer  species:  75  species  and  529 
reactions  versus  155  species  and  689  reactions  in  the  Marinov  mechanism.  Fewer  species 
and  reactions  mean  that  the  reduced  mechanism  is  likely  to  run  faster  and  to  exhibit  better 
numerical  behavior  and  agreement  with  the  detailed  mechanism  using  fewer  non-QSS 
species.  Thus,  the  Wang  mechanism  was  selected  as  the  basis  for  the  ethylene  reduced 
mechanisms  to  be  used  in  this  work. 

4.2  Complex  Fuels 
4.2.1  Heptane 

Several  chemical  kinetic  mechanisms  for  n-heptane  combustion  exist  in  the  literature. 

One  commonly  used  « -heptane  mechanism  is  that  of  Curran  et  al.  (1998),  which  is 
intended  to  cover  the  entire  range  of  conditions  from  low-temperature  (600  to  900  K) 
pyrolysis  and  oxidation  to  high-temperature  combustion.  During  Phase  I,  a  subset  of  this 
mechanism  containing  only  those  reaction  steps  and  species  that  are  only  important  at 
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higher  temperatures  (>900  K)  was  used  as  a  basis  for  reduction.  This  high-temperature 
subset  is  still  a  large  mechanism,  with  105  species  and  808  elementary  steps. 

In  contrast,  the  n-heptane  mechanism  of  Held  et  al.  (1997)  is  considerably  simpler.  This 
mechanism  compares  well  to  experiments  with  comparatively  very  few  species  by 
empirically  modeling  the  initial  fuel  breakdown.  The  Curran  mechanism  contains 
reaction  steps  for  abstraction  of  H  atoms  from  the  fuel  by  a  number  of  radical  species  to 
form  several  heptyl  radical  isomers.  These  heptyl  radicals  undergo  isomerization 
reactions  and  decompose  through  a  number  of  routes  to  form  various  species  with  two  to 
five  carbon  atoms.  Held  et  al.  achieve  considerable  simplification  by  bypassing  the 
formation,  isomerization,  and  breakdown  of  the  various  heptyl  radicals  and  by  allowing 
the  n-heptane  fuel  to  decompose  directly  into  smaller  reaction  products,  often  with  three 
or  four  products  on  the  right-hand  side  of  an  elementary  reaction  step.  Because  of  its 
good  comparison  to  experiments  and  comparative  simplicity,  the  Held  mechanism  was 
chosen  for  reduction. 

4.2.2  JP-8 

The  primary  candidates  for  a  detailed  mechanism  for  JP-8  combustion  are  the 
mechanisms  of  Violi  et  al.  (2002)  and  Mawid  and  Sekar  (2002).  We  have  chosen  to  use 
the  Violi  mechanism  because  of  our  access  to  it  and  because  it  compares  very  well  to  the 
limited  data  available  on  JP-8  and  kerosene  combustion.  This  mechanism  was  created  by 
collaboration  between  the  University  of  Utah  and  Italian  researchers.  One  of  the  authors 
of  the  mechanism,  Prof.  Adel  Sarofim,  is  a  technical  advisor  to  REI  on  this  project, 
allowing  us  to  actively  engage  in  discussions  regarding  the  development  and  use  of  this 
mechanism.  The  Violi  mechanism  was  not  originally  in  the  Chemkin  format  needed  for 
reduction  by  CARM.  The  Violi  mechanism  was  converted  to  Chemkin  format  by  REI  as 
part  of  this  project. 

Experimental  results  (Doute  et  al.,  1995)  for  kerosene  combustion  in  a  fuel-rich 
(equivalence  ratio=1.7)  PFR  have  been  modeled  successfully  using  the  Violi  mechanism. 
The  comparison  of  the  measured  and  calculated  profdes  is  shown  in  Figure  3,  where  the 
surrogate  is  composed  of  73.5  mol%  n-dodecane,  5.5  mol%  iso-octane,  10  mol% 
methylcyclohexane,  1 1  mol%  toluene,  and  1  mol%  benzene.  The  mechanism  reproduces 
the  profdes  of  fuel  components,  products,  and  many  intermediate  species  quite  well. 

Figure  4  compares  ignition  delay  predictions  made  by  the  Violi,  and  Mawid  and  Sekar 
mechanisms  with  measurements  from  Mullins  (1955)  and  Freeman  and  Lefebvre  (1984). 
The  Violi  mechanism  uses  a  six-specie  surrogate  blend  to  represent  a  typical  batch  of  JP- 
8  fuel.  The  Mawid  and  Sekar  mechanism  uses  12  species.  Having  more  species  in  the 
surrogate  is  a  disadvantage  for  mechanism  reduction  because  each  species  in  the 
surrogate  and  some  intermediates  related  to  it  must  be  kept  in  the  reduced  mechanism. 
Though  the  Violi  mechanism  has  not  been  tuned  to  match  these  data,  the  agreement  is 
somewhat  better  than  that  of  Mawid  and  Sekar  over  the  limited  data  range. 
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ignition  delay  (sec) 


Figure  2.  Comparison  of  Ignition  Delay  as  a  Function  of  Temperature  for 
Experiments,  Detailed  Chemistry  and  the  Reduced  Mechanism  for  Ethylene-air 
Mixtures  at  E.R.s  of  0.5, 1.0,  and  2.0. 
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Figure  3.  Comparison  of  the  Violi  Mechanism’s  Predictions  to  Experimental 
Results  in  a  Kerosene-air  PFR.  T=1160  K,  P=1  atm.,  E.R.=1.7. 
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Figure  4.  Comparison  of  the  Mechanism  Predictions  for  JP-8  (Freeman  and 
Lefebvre,  1984)  or  Kerosene  (Mullins,  1955)  Ignition  Delay.  P  =  1  atm,  E.R.  =  0.5. 
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5.0  Development  of  Reduced  Mechanisms 

In  this  section,  the  reduced  mechanisms  that  have  been  created  using  CARM  are 
described  and  compared  to  detailed  chemistry  and,  where  available,  experimental  data. 

Of  primary  interest  for  scramjet  combustor  simulations  are  the  ignition  delay  and  heat 
release  rate.  Ignition  delay  is  very  sensitive  to  the  details  of  the  chemical  kinetics,  while 
the  overall  heat  release  is  a  function  of  the  thermodynamic  properties  of  the  fuel  and 
major  products.  Ignition  delay  for  reduced  and  detailed  chemistry  has  been  compared 
over  a  range  of  stoichiometries  and  temperatures  for  all  the  reduced  mechanisms  created 
and  implemented  in  this  project.  It  is  believed  that  if  the  ignition  delay  is  correctly 
modeled,  then  the  heat  release  rates  and  distributions  of  important  species  will  be 
modeled  accurately  as  well. 

5.1  CARM  Phase  I  Background 

There  are  four  basic  steps  in  the  formulation  of  a  reduced  chemical  kinetic  mechanism:  1) 
identification  of  the  appropriate  detailed  mechanism  containing  the  essential  species  and 
elementary  reaction  steps,  2)  identification  of  appropriate  quasi-steady-state  (QSS) 
approximations,  3)  elimination  of  reactions  through  use  of  the  algebraic  relations 
obtained  in  step  2,  and  4)  solution  of  the  coupled  and  nonlinear  set  of  algebraic  equations 
obtained  in  the  previous  steps  to  find  the  reaction  rates  of  the  remaining  species.  CARM 
automates  steps  2  through  4,  producing  source  code  for  the  calculation  of  the  chemical 
source  terms  defined  by  the  reduced  mechanism.  As  inputs,  CARM  uses  a  set  of  test 
problem  results  representing  conditions  of  interest  to  rank  species  by  the  error  introduced 
by  assuming  they  are  in  steady  state.  The  subroutine  produced  by  CARM  contains  code 
that  iteratively  solves  the  coupled,  nonlinear  set  of  algebraic  equations  giving  the 
concentrations  of  the  quasi-steady-state  species.  These  concentrations  are  used  along 
with  the  rates  of  the  elementary  reactions  from  the  detailed  mechanism  to  calculate  the 
chemical  source  terms  for  the  non-QSS  species. 

During  Phase  I,  the  ability  of  the  CARM  software  to  create  reduced  mechanisms  for 
hydrocarbon  fuels  (ethylene  and  « -heptane  using  two  different  mechanisms)  was 
demonstrated.  The  ability  of  reduced  mechanisms  created  using  CARM  to  correctly 
predict  temperatures  and  species  concentrations  in  simple  reactors  calculations  (Plug 
Flow  Reactors  (PFRs)  and  Perfectly  Stirred  Reactors  (PSRs))  was  examined.  Before  the 
Phase  I  work,  CARM  had  only  been  demonstrated  on  methane. 

5.2  CARM  Software  Improvements 

During  this  project,  several  improvements  were  made  to  CARM.  These  improvements 
have  been  aimed  at  improving  the  speed  and  robustness  of  the  subroutines  produced  by 
CARM,  as  well  as  allowing  for  nonstandard  features  of  some  detailed  mechanisms  of 
interest,  such  as  the  Violi  mechanism. 

5.2.1  Chemical  Rate  Table 

Chemical  rate  expressions  contain  exponential  functions  of  temperature  only.  Evaluation 
of  exponential  functions  can  be  costly  for  repeated  calculations  as  required  by  most  CFD 
calculations.  In  the  latest  version  of  CARM,  a  new  option  is  included  for  building  a 
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chemical  rate  table  as  a  function  of  temperature.  This  lookup  table  is  built  once.  The 
chemical  rates  are  determined  from  an  efficient  linear  interpolation  scheme.  Nearly  a 
factor  of  three  speedup  in  the  rate  expression  computations  has  been  observed. 

5.2.2  Acceptance  of  Nonstandard  Mechanism  Formatting 

CARM  was  modified  to  work  for  detailed  mechanisms  that  do  not  conform  to  standard 
Chemkin  format.  CARM  can  now  operate  on  detailed  mechanisms  that  have  noninteger 
stoichiometric  coefficients  and  up  to  12  products  on  the  right-hand  side  of  an  elementary 
reaction  step.  This  modification  was  necessary  for  reduction  of  the  Violi  et  al.  (2002)  JP- 
8  mechanism. 

5.2.3  FORTRAN  90  Subroutine  Output 

CARM  has  been  updated  to  provide  a  FORTRAN-90  compatible  reduced  mechanism 
source  subroutine.  This  capability  has  been  extensively  tested.  FORTRAN-90 
subroutines  are  significantly  easier  to  implement  into  VULCAN. 

5.2.4  Newton  Iteration  for  QSS  Species 

Reduced  mechanisms  produced  by  previous  CARM  versions  used  a  fixed-point  iteration 
scheme  to  solve  for  the  QSS  species  concentrations.  Occasional  convergence  failures 
have  been  observed  when  a  limit  is  placed  on  the  number  of  fixed-point  iterations.  This 
failure  sometimes  slows  computation,  and  can  impair  the  accuracy  of  the  overall 
calculation.  Since  the  fixed-point  iteration  scheme  has  a  slow  convergence  rate  on  the 
order  of  unity  (i.e.,  almost  linear),  modifications  to  CARM  have  been  made  to  provide 
the  capability  of  performing  Newton  iterations  when  the  fixed-point  scheme  fails.  The 
major  coding  effort  of  these  modifications  focused  on  automatic  output  of  FORTRAN 
code  to  evaluate  the  Jacobian  terms  needed  in  the  Newton  iteration.  Matrix  inversion  is 
the  most  expensive  part  of  Newton  iteration;  a  Lower-Upper  (LU)  matrix  decomposition 
scheme  is  used.  The  automatically  generated  code  is  verified  by  comparisons  of  results 
obtained  from  a  PSR  model. 

Several  reduced  mechanisms  with  different  numbers  of  steady-state  species  have  been 
developed  for  assessment  of  the  performance  of  the  combined  scheme.  Table  1  compares 
the  ratio  of  CPU  times  used  by  Newton  and  fixed-point  iteration  along  with  the 
percentage  of  CPU  time  used  for  matrix  inversion.  As  expected,  the  CPU  time  required 
for  Newton  iteration  is  larger  than  that  required  by  the  fixed-point  scheme  and  the 
majority  of  the  time  is  used  for  matrix  inversion.  Matrix  inversion  uses  over  90  percent  of 
the  CPU  time  when  the  total  number  of  steady-state  species  exceeds  20. 


15 


Table  1.  Comparison  of  Ratio  of  CPU  Times  Used  by  Newton  Iteration  and  by  the 
Fixed-point  Iteration  Along  with  the  Percentage  of  CPU  Time  Used  for  Matrix 
Inversion. 


Reduced 
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Number  of  steady- 

CPU  Newton/ 

Percent  CPU  for 

mech. 

of  species 

state  species 

CPU  fixed-point 

matrix  inversion 

H2+  NOx 

11 

11 

15.9 

61.6 

Skeletal  CH4 

5 

11 

6.79 

83.8 

GRI  2. 11 

12 

33 

37.2 

92.2 

GRI  3.0 

42 

10 

3.4 

84.6 

GRI  3.0 

32 

20 

19.2 

92.5 

GRI  3.0 

19 

33 

27.9 

94.5 

GRI  3.0 

12 

40 

33.9 

92.1 

Figure  5  presents  the  CPU  times  used  by  the  matrix  inversion  versus  the  number  of 
variables  (the  steady-state  species)  showing  the  expected  cubic  dependence.  Figures  6 
and  7  plot  the  ratio  of  CPU  times  used  by  Newton  and  by  fixed-point  iteration  based  on  a 
PSR  calculation  on  linear  and  log-log  scales,  respectively.  Due  to  the  CPU  times  used  by 
other  parts  of  the  PSR  code,  the  dependence  of  total  CPU  time  on  the  number  of  steady- 
state  species  scales  roughly  as  ~(NSS)1<5,  where  Nss  stands  for  the  total  number  of  steady- 
state  species.  Figure  8  compares  the  number  of  iterations  required  to  reach  convergence 
by  the  two  schemes.  As  can  be  seen  in  the  figure,  the  fixed-point  scheme  requires  a  large 
number  of  iterations  compared  to  the  Newton  method.  On  average,  the  fixed-point 
scheme  performs  22.6  iterations  to  reach  convergence,  while  the  Newton  scheme  requires 
only  3.01  iterations.  However,  due  to  the  large  CPU  time  required  by  the  Newton 
method,  the  overall  CPU  time  is  much  larger  than  that  used  by  the  fixed-point  scheme. 
This  conclusion  is  somewhat  expected;  the  main  motivation  for  performing  Newton 
iteration  is  in  case  of  convergence  failure  by  the  fixed-point  scheme. 

Application  of  combined  fixed  point  and  Newton  iterations  for  steady  state  species  was 
carried  out  using  a  well  mixed  reactor  model.  The  reduced  chemistry  is  based  on  Gas 
Institute  Research  (GRI)-3.0  (GRI  Mech  web  site)  CH4-air  combustion  with  15  steps  and 
33  steady-state  species.  The  stiff  ordinary  differential  equation  (ODE)  integrator 
DDASAC  (Double  precision  Differential  Algebraic  Sensitivity  Analysis  Code)  made  a 
total  of  38,577  calls  to  the  subroutine  CKWYP;  2.2  percent  of  these  calls  with  the  fixed- 
point  iterations  failed  to  converge  within  20  iterations.  The  average  number  of  fixed  point 
iteration  was  6.3  while  the  average  number  of  Newton  iteration  was  3.2.  Figure  9 
illustrates  a  sample  set  of  calling  statistics  showing  that  when  fixed-point  iteration 
exceeds  the  assigned  maximum  limit,  Newton  iteration  is  called,  leading  to  fast 
convergence.  Although  Newton  iteration  is  expensive,  many  of  the  convergence  failures 
of  the  fixed-point  iteration  method  can  be  avoided. 

An  exploration  run  was  conducted  with  the  fixed-point  iterations  only.  During  this  run, 
the  stiff  ODE  integrator  DDASAC  made  46,143  calls  to  the  subroutine  CKWYP,  about 
20%  more  than  using  the  combined  scheme.  This  is  due  to  the  fact  that  more  calls  are 
needed  when  iteration  to  calculate  the  steady-state  species  fails  to  converge.  It  is 
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interesting  that  the  actual  CPU  runtime  with  the  fixed  iteration  scheme  is  about  1.3% 
more  expensive  than  the  combined  scheme.  This  may  be  fortuitous  under  this  special 
application  as  Newton  iteration  is  more  expensive.  However,  the  speeding  convergence 
offered  by  the  Newton  iteration  can  offset  the  computer  time  by  having  fewer  overall 
iterations.  This  development  provides  a  robust  iteration  scheme  to  achieve  a  converged 
solution  when  the  fixed-point  iteration  scheme  fails. 


Figure  5.  CPU  Times  Used  in  Inverting  a  Square  Matrix  Versus  Number  of  Steady- 
State  Species  Showing  the  Statistics  with  the  Expected  Cubic  Dependence. 
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Figure  6.  Ratio  of  Total  CPU  Time  versus  Number  of  Steady-State  Species  on  a 
Linear  Scale. 
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Figure  7.  Ratio  of  Total  CPU  Time  versus  Number  of  Steady-State  Species  on  a 
Log-log  Scale. 


Figure  8.  Comparison  of  the  Number  of  Iterations  to  Reach  Converged  Solutions 
for  the  Fixed-Point  Iteration  Scheme  and  the  Newton  Iteration  Method  for  a  PSR 
Calculation. 
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Figure  9.  Number  of  Iterations  in  a  Combined  Fixed-Point  and  Newton  Iteration 
Scheme. 

5.2.5  Skeletal  Mechanism  Development 

A  skeletal  mechanism  is  a  smaller  version  of  a  detailed  chemical  kinetic  mechanism  with 
unimportant  species  and  reactions  deleted.  For  fuels  requiring  very  large  mechanisms, 
creation  of  a  good  skeletal  mechanism  is  an  important  step  in  creating  the  reduced 
mechanism  with  CARM.  Approaches  for  constructing  accurate  skeletal  mechanisms 
were  evaluated  and  new  improvements  to  CARM  were  developed.  Three  different 
methods  including  detailed  reduction  (DR)  (Wang  and  Frenklach,  1991),  atomic  flux 
analysis  (Soyhan  et  al.,  2001),  and  the  original  CARM  scheme  (Chen,  1997),  were 
assessed  for  development  of  accurate  skeletal  mechanisms,  especially  for  predictions  of 
autoignition  delay.  Two  features  in  the  original  CARM  method  were  identified  to  be 
undesirable.  One  is  the  arbitrary  threshold  used  for  the  initial  selection  of  species  to  be 
kept  in  the  skeletal  mechanism.  The  other  one  is  the  required  sensitivity  information  that 
is  known  to  be  time  consuming  to  obtain,  especially  for  large  mechanisms.  Furthermore, 
the  large  memory  needed  for  storing  sensitivity  information  presents  a  limiting  factor  on 
the  number  of  combustion  states  that  can  be  included  during  the  development.  Efforts 
have  been  made  to  remove  these  two  undesirable  features  in  the  current  CARM  version 
by  considering  alternatives  (Tham  and  Chen,  2003). 

The  atomic  flux  analysis  method  was  explored  as  an  alternative  to  the  previously  existing 
CARM  approach  for  selecting  the  species  to  be  kept  in  the  skeletal  mechanism.  However, 
experience  with  large  mechanisms,  such  as  iso-octane,  revealed  that  the  atomic  flux 
analysis  can  be  quite  time  consuming  and  is  unsuitable  for  running  CARM  interactively. 
Another  alternative  is  reaction  matrix  analysis  (Bendtsen,  et  al.,  2001),  which  is  similar  to 
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atomic  flux  analysis  but  requires  far  less  computation  time  for  large  mechanisms.  Instead 
of  a  user-specified  threshold  value  to  distinguish  species  importance,  a  new  method  of 
tracing  the  most  important  reaction  paths  of  major  reactants  and  products  is  employed  in 
the  reaction  matrix  analysis.  Our  study  of  the  reaction  matrix  analysis  technique  reveals 
that  similar  results  to  those  from  atomic  flux  analysis  are  obtained  with  much  less 
computation  time. 

The  DR  approach  focuses  on  the  reaction  rate  and  heat  release  rate  of  each  step.  On  the 
other  hand,  the  original  CARM  approach  focuses  on  sensitivity  information  and 
construction/destruction  rate  of  each  species.  The  new  version  of  CARM  includes 
advantageous  features  from  different  approaches.  The  hybrid  approach  retains  the  heat 
release  rate  analysis  from  the  DR  approach  along  with  the  construction/destruction  rate 
analysis  from  the  CARM  approach  to  efficiently  produce  a  more  accurate  skeletal 
mechanism.  The  initial  selection  of  species  is  performed  using  the  newly  developed 
tracing  method. 

Assessment  of  the  performance  of  the  latest  CARM  version  is  conducted  for  an 
autoignition  study  with  a  detailed  hydrogen  mechanism.  Three  methods,  including  DR, 
the  original  CARM,  and  the  newly  modified  CARM,  are  used  to  create  skeletal 
mechanisms  with  the  constraint  that  they  contain  the  same  number  of  species. 
Comparisons  of  predicted  temperature  and  H2O2  concentration  during  autoignition  are 
presented  in  Figure  10.  As  seen  in  the  comparisons,  the  original  CARM  approach 
predicts  lower  temperatures  near  the  end  of  autoignition.  The  DR  approach  miscalculates 
the  H2O2  mole  fraction.  The  newly  improved  CARM  approach  (hybrid  approach)  gives  a 
skeletal  mechanism  that  agrees  closely  with  the  detailed  mechanism  in  terms  of  both 
temperature  and  H2O2.  This  new  capability  has  proven  to  be  useful  for  creating  reduced 
mechanisms  for  large  hydrocarbons  and  blends  representing  real  fuels. 


Figure  10.  Comparisons  of  Three  Different  Approaches  for  Creating  Skeletal 
Mechanisms. 
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5.3  Reduced  Mechanisms  for  Simple  Fuels 

5.3.1  Reduced  Mechanism  for  Hydrogen-Air  Combustion 

Using  CARM,  two  reduced  mechanisms  designed  to  model  ignition  of  hydrogen-air 
mixtures  were  created.  The  reduced  mechanisms  were  based  on  H2-O2  reactions 
extracted  from  a  large  hydrocarbon  mechanism  (Curran  et  ah,  1998).  The  reduced 
mechanisms  reduce  the  number  of  species  from  nine  in  the  detailed  mechanism  to  six  or 
seven,  by  assuming  that  OH  and  H2O2,  or  OH,  H2O2,  and  HO2  are  in  QSS.  These  reduced 
mechanisms  have  been  tested  against  detailed  chemistry  in  constant  pressure  ignition 
delay  calculations  for  pressure  =  1  atm,  E.R.  =  0.5,  1.0,  and  2.0,  and  initial  temperature 
=  1000  to  1400  K.  Some  of  the  results  are  plotted  in  Figure  11.  The  seven-specie  reduced 
mechanism  works  very  well.  The  six-specie  mechanism  results  in  some  error  in  the 
computed  ignition  delay  at  the  lower  end  of  the  temperature  range  QSS  for  OH  and  H2O2 
makes  almost  no  difference,  but  error  begins  to  appear  when  HO2  is  assumed  in  QSS. 

The  seven-species  reduced  mechanism  shows  excellent  agreement  with  detailed 
chemistry  and  was  selected  for  initial  implementation  into  VULCAN.  Reduction  from 
nine  to  seven  species  is  not  a  large  reduction,  but  allowed  the  VULCAN  implementation 
to  be  tested  on  a  relatively  simple  case  for  which  experimental  data  exist. 


Figure  11.  Comparison  of  Ignition  Delay  as  a  Function  of  Temperature  for  Detailed 
and  Reduced  Mechanisms  for  H2-Air  Mixture  at  E.R.s  of  0.5  and  1.0,  Respectively, 

P  =  1  atm. 

5.3.2  Reduced  Mechanisms  for  Ethylene-Air  Combustion 

Ethylene  (C2H4),  is  a  small  hydrocarbon  which  has  been  used  previously  in  combustor 
simulations  and  experiments  as  a  scramjet  fuel  (Baurle  and  Eklund,  2002).  Based  on  the 
ethylene  detailed  mechanism  comparisons,  reduced  mechanisms  were  created  and  tested 
based  on  the  Wang  et  al.  (1999)  mechanism.  A  reduced  mechanism  with  20  species  (cf. 
75  in  the  detailed  mechanism)  gives  excellent  agreement  with  detailed  chemistry  and 
experimental  results  (see  Figure  12).  A  15-specie  reduced  mechanism,  also  based  on  the 
Wang  detailed  chemistry,  did  not  perform  satisfactorily,  giving  ignition  delay  errors  for 
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some  conditions  of  an  order  of  magnitude  or  more.  The  20-specie  reduced  mechanism 
has  been  implemented  into  VULCAN  to  further  evaluate  the  numerical  accuracy  and 
robustness  of  the  implementation. 

The  non-QSS  species  in  this  mechanism  are  H2,  H,  O,  O2,  OH,  H2O,  HO2,  CH3,  CH4,  CO, 
C02,  HCO,  CH20,  C2H2,  C2H3,  C2H4,  HCCO,  CH2CO,  CH3CHO,  and  N2.  This  new 
reduced  mechanism  is  an  improvement  on  previously  published  reduced  mechanisms  for 
ethylene  combustion  (Montgomery  et  ah,  2002b),  being  based  on  a  more  up-to-date 
detailed  mechanism  and  tuned  for  ignition  delay  under  conditions  of  interest  for  a 
scramjet  combustor. 


Figure  12.  Comparison  of  Ignition  Delay  as  a  Function  of  Temperature  for 
Experiments,  Detailed  Chemistry  and  the  Reduced  Mechanism  for  Ethylene-Air 
Mixtures. 

5.4  Reduced  Mechanisms  for  Complex  Fuels 

5.4.1  Reduced  Mechanism  for  A-Heptane-Air  Combustion 

A  19-specie  n-heptane  reduced  mechanism  based  on  the  detailed  mechanism  of  Held  et 
al.  (1997)  was  created  using  CARM.  Figure  13  compares  calculated  ignition  delays  for  a 
range  of  E.R.s  and  temperatures  for  detailed  and  reduced  chemistry.  Despite  the  fact  that 
/? -heptane  is  a  larger  fuel  with  a  more  complex  combustion  process  than  ethylene,  a 
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smaller  reduced  mechanism  (19  species  cf.  20  for  ethylene)  could  be  created  with  similar 
accuracy  due  to  the  reductions  inherent  in  the  semi-empirical  mechanism  of  Held  et  al. 
(1997).  The  species  contained  in  the  reduced  mechanism  are  n- C7H16,  02,  C2H4,  CO, 
C02,  C3H6,  CH4,  C2H6,  C4H8,  C5H10,  C6H12,  C2H2,  C3H5,  CH20,  H,  H202,  H2,  H20,  and 

n2. 


1/T  (1/K) 


1/T  (1/K)  1/T  (1/K) 

Figure  13.  Ignition  Delay  Calculations  for  n-Heptane  in  Air  for  Detailed  Chemistry 
and  the  19-Specie  Reduced  Mechanism. 

5.4.2  Reduced  Mechanisms  for  JP-8-Air  Combustion 

JP-8  fuel  is  a  complex  mixture  of  thousands  of  hydrocarbon  components  that  is 
manufactured  to  meet  certain  property  standards.  The  exact  composition  varies  in  time 
and  place  depending  on  crude  feedstocks  and  refinery  processes.  Modeling  JP-8 
combustion  requires  selection  of  a  surrogate  fuel  blend  containing  a  few  species  with 
known  kinetics  that  reproduces  the  essential  combustion  behavior  of  a  typical  batch  of 
the  actual  fuel  (Edwards  &  Maurice,  2001).  The  newly  upgraded  CARM  (with 
modifications  described  in  Section  5.2.2  )  was  used  for  development  of  a  JP-8  reduced 
mechanism  based  on  the  semi-empirical  mechanism  developed  by  Violi  et  al.  (2002) 
using  a  six-species  surrogate  JP-8  fuel  mixture,  which  is  listed  in  Table  2.  This  surrogate 
blend  was  developed  to  match  the  boiling  curve  of  a  batch  of  JP-8  and  not  for  any 
combustion  characteristics,  yet  was  shown  earlier  to  give  excellent  agreement  to  ignition 
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delay  data  (see  Figure  4).  The  semi-empirical  mechanism  contains  216  species  and  4,826 
steps. 


Table  2.  JP-8  Surrogate  Blend  (Violi  et  al.  2002). 


Compound 

Formula 

Mole  Fraction 

dodecane 

C12H26 

0.30 

tetradecane 

Ci4H3o 

0.20 

iso-octane 

c8h10 

0.10 

methyl-cylcohexane 

c7h14 

0.20 

m-xylene 

c8h10 

0.15 

tetralin 

C10H12 

0.05 

A  transient  PSR  was  run  with  initial  temperature  of  1200  K  and  pressure  of  1  atm  to 
provide  the  needed  information  for  development  of  skeletal  mechanisms.  Several  trial 
versions  with  different  sizes  were  tested,  and  the  final  version  contains  147  species  and 
833  steps.  Further  reduction  using  CARM  was  carried  out  to  develop  reduced  mechanism 
of  various  sizes.  As  the  choice  of  the  QSS  species  pool  can  impact  the  stiffness  of  the 
reduced  chemistry,  CARM  was  run  iteratively  with  manageable  QSS  pools  to  meet  the 
dual  criteria  of  high  accuracy  without  unmanageable  stiffness.  A  35-specie  reduced 
mechanism  was  found  to  be  quite  accurate  for  a  JP-8  surrogate  fuel  with  six  components. 
The  non-QSS  species  contained  in  the  35-specie  mechanism  are  Fh,  CH4,  C2H2,  m- 
xylene,  C2H4,  indenyl,  C2H6,  tetralin,  H,  O,  p- C3H4,  CH3,  a-C 3H4,  C3H6,  C12H25, 
butadiene,  butene,  isobutane,  C5H10,  methyl-cylcohexane,  H2O,  CO,  CO2,  O2,  OH,  H2O2, 
C7H14,  CH2O,  CH2CO,  benzaldehyde,  iso-octane,  dodecane,  tetradecane,  C8H6,  and  N2. 

A  33-specie  reduced  mechanism  did  not  give  sufficient  accuracy.  Figure  14  presents 
comparisons  of  the  predicted  ignition  delay  times  among  the  detailed  and  reduced 
mechanisms  showing  excellent  agreement.  This  reduced  mechanism  performs 
significantly  better  than  those  created  previously  (Montgomery  et  al.,  2002c)  from  the 
Mawid  and  Sekar  mechanism. 
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Figure  14.  Comparison  of  Calculated  Ignition  Delays  for  JP-8  in  Air. 

In  an  attempt  to  create  a  reduced  mechanism  for  JP-8  having  significantly  fewer  than  35 
species,  further  simplification  was  attempted.  The  Violi  et  al.  (2002)  mechanism  uses  a 
six-species  fuel  surrogate  blend  to  represent  JP-8.  It  was  hypothesized  that  a  significantly 
smaller  reduced  mechanism  could  be  created  by  using  a  simpler  fuel  surrogate.  Each 
species  in  the  fuel  surrogate  plus  some  of  its  breakdown  products  must  be  included  in  the 
reduced  mechanism.  Through  a  trial-and-error  process,  we  have  created  a  two-species 
surrogate  containing  the  normal  dodecane  ( /1-C12H26 )  and  m-x ylene  (CsHio)  that  gives 
calculated  ignition  delay  results  that  are  indistinguishable  from  those  found  using  the  six- 
species  surrogate  as  shown  in  Figure  15.  The  two  species  in  the  surrogate  represent, 
respectively,  the  paraffin  and  aromatic  components  of  the  actual  fuel. 
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Figure  15.  Comparison  of  Calculated  Ignition  Delays  Using  Six-  and  Two-Species 
Fuel  Surrogates  for  JP-8. 

However,  satisfactory  reduced  mechanisms  for  the  2-specie  surrogate  having  no  fewer 
than  3 1  species  could  be  created,  not  a  particularly  significant  savings  in  comparison  to 
the  35-specie  mechanism.  Therefore,  it  was  decided  to  proceed  with  the  implementation 
of  the  35-specie  reduced  mechanism  into  the  2-D  flame  case  in  VULCAN.  The  35- 
specie  mechanism  has  been  verified  over  a  wide  temperature  range  and  for  E.R.s  ranging 
from  nearly  pure  air  to  nearly  pure  fuel  in  a  PFR.  The  mechanism  does  not  present  any 
numerical  difficulties  under  these  conditions.  This  reduced  mechanism  and  surrogate 
blend  give  an  excellent  representation  of  the  components  of  actual  JP-8  fuel  and  their 
important  combustion  characteristics,  such  as  heat  release  rate  and  ignition  delay. 
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6.0  In  Situ  Adaptive  Tabulation  (ISAT) 

The  CARM  produced  reduced  mechanism  may  require  a  great  deal  of  CPU  time  to 
perform  the  iterative  calculations  to  obtain  the  species  reaction  rates,  depending  on  the 
number  of  steady-state  species.  For  example,  in  a  CFD  simulation  of  a  supersonic  C2H4 
diffusion  flame  with  a  20-specie  reduced  mechanism  (to  be  detailed  later),  more  than 
98%  of  the  total  CPU  time  is  consumed  by  the  reduced  mechanism  subroutine  if  the 
reaction  rates  are  evaluated  directly.  In  order  to  account  for  turbulence-chemistry 
interactions,  evaluation  of  the  turbulent  reaction  rates  involves  numerous  calls  to  the 
reduced  mechanism  subroutine,  which  is  even  more  expensive.  Thus,  significant  CPU 
savings  can  be  achieved  by  tabulating  the  chemical  source  terms. 

Direct  tabulation  is  efficient  and  accurate  when  there  are  only  2  or  3  independent 
variables.  However,  it  is  not  feasible  for  a  reduced  mechanism  that  contains  10  or  more 
species.  In  situ  adaptive  tabulation  (ISAT),  developed  in  the  context  of  probability 
density  function  (PDF)  methods  (Pope,  1997),  is  reported  to  be  able  to  achieve  a  speed¬ 
up  factor  of  1000. 

The  “accessed  region”  of  the  composition  space  is  defined  as  the  set  of  all  compositions 
that  occur  in  a  reactive  flow  (or  calculation  of  this  flow).  The  basic  idea  of  ISAT  is  based 
on  a  crucial  observation:  the  accessed  region  is  much  smaller  than  the  realizable  region.  It 
is  sufficient  to  tabulate  only  the  accessed  region  rather  than  the  entire  realizable  region. 
An  unstructured  adaptive  table  is  built  up  during  the  reactive  flow  simulation,  which  is 
referred  as  in  situ  tabulation. 

In  this  research  program,  a  scheme  developed  for  the  PDF  method  (Chen  et  al.,  2000) 
was  modified  and  implemented  into  VULCAN.  The  details  of  the  development  and 
implementation  are  discussed  in  the  next  two  sections. 

6.1  Mathematical  Formulation 

The  ISAT  algorithm  employed  here  follows  closely  the  principles  proposed  by  Pope 
(1997),  except  that  the  target  functions  considered  here  are  the  species  reaction  rates 
instead  of  the  species  concentrations  after  a  given  time  step.  The  reaction  rates  can  be 
written  as  functions  of  mass  fractions  of  species,  temperature  and  pressure: 

,  t(Yx,Y2,...,Ym,T ,  p)  i  l,...,  ns,  (1) 

where  ns  is  the  number  of  species.  More  generally,  in  a  vector  form,  we  have  R=R(Z) 
where  R=(Ri,...,Rn)t=  (  1,...,  ns)T  and  Z=(Zi,...,ZM)T=(Yi,...,Yns,T,p)T 

The  interpolation  formula  used  in  ISAT  is 

R(Z°  Z)  R(Z°)  ^  Z  ,  (2) 

where  Z°is  the  tabulation  point,  and  (Z°  Z)  is  the  query  point. 
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(3) 


The  components  of  the  Jacobian  matrix,  r(z0)/  z,  are  calculated  numerically: 

^(Z,0,...^)  -  R(Z«,...^  Rtf,. 

i  1  ,..N,  j  1  ,...M 

Note  that  evaluation  of  the  Jacobian  matrix  is  expensive.  It  requires  (M+l)  times  as  much 
CPU  as  the  evaluation  of  the  reaction  rate  vector. 

The  interpolation  error  in  ISAT  is  directly  associated  with  the  magnitude  of 
R  (  Rm)t  .  The  error  is  controlled  by  setting  a  tolerance  for  all  the  components 
of  R  .  The  ellipsoid  of  accuracy  (EOA)  of  a  tabulation  point  z°  is  defined  as  the  domain 
surrounding  the  point  such  that 

I  J?  .  I  , . , 

Max  —J—  tol  j  1  ,...,ns,  (4) 

jfVy',nom 

where  Rj  nom  is  the  value  for  normalization,  and  tol  is  the  tolerance. 

6.2  Dynamic  Tabulation 

The  table  for  ISAT  is  built  as  the  calculation  proceeds.  Initially  the  table  is  empty.  On 
receipt  of  the  first  query  Z? ,  the  first  table  entry  is  generated.  As  it  receives  the  second 
query,  the  ISAT  subroutine  checks  first  if  the  query  lies  within  the  EOA  of  the  existing 
table  entry  for  an  interpolation.  If  it  is  not,  the  second  entry  is  generated.  A  cutting  plane 
is  also  generated  to  divide  the  composition  space  into  two  parts.  In  our  study,  the  cutting 
plane  is  defined  as  the  plane  perpendicular  to  the  vector  connecting  the  two  tabulation 
points,  as  given  in  the  following: 

Z  ^(Z1  Z2)  (Zj  Z2)  0  ,  (5) 

where  Zi  and  Z2  are  the  two  tabulation  points.  A  sketch  of  the  cutting  plane  in  a  2D  space 
is  shown  in  Figure  16. 


Figure  16.  Sketch  of  Cutting  Plane  in  Relation  to  the  First  and  Second  ISAT  Table 
Entries. 
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For  the  third  query  point,  it  is  first  decided  which  side  of  the  cutting  plane  the  query  point 
belongs  to.  Then  the  point  is  checked  to  determine  whether  the  query  is  close  enough  to 
the  table  entry  on  the  same  side.  If  it  is,  interpolation  is  performed;  if  it  is  not,  a  new  entry 
is  generated  and  a  new  cutting  plane  is  defined  to  separate  the  old  entry  and  the  new 
entry.  The  same  algorithm  repeats  as  more  query  points  are  received. 


The  dynamic  table  is  stored  in  a  binary  tree,  as  sketched  in  Figure  17.  Each  leaf  of  the 
tree  represents  an  entry  of  the  table,  where  the  tabulation  point  Z° ,  the  reaction  rates 

R(Z°) 


R(Z  )  ,  and  the  Jacobian  matrix 


are  stored.  Each  node  of  the  tree  represents  a 


cutting  plane. 


Figure  17.  Sketch  of  the  Binary  Tree.  Each  leaf  (black  dot)  contains  a  table  entry; 
each  node  (circle)  contains  information  about  the  cutting  plane. 

Each  time  a  query  point  is  received,  the  binary  tree  is  traversed  until  a  leaf  is  reached. 
Ideally,  one  would  like  to  find  the  leaf  Z°  that  is  closest  to  query  point  Z?  in  the  sense 
that  the  linear  approximation  error  is  minimized.  However,  it  is  computationally 
expensive  to  find  the  closest  leaf.  The  binary  search  used  instead  is  computationally 
inexpensive,  and  yields  a  leaf  Z°  that  is  likely  to  be  close  to  Z1' . 
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7.0  Turbulence- Chemistry  Interaction 

The  modeling  of  the  mean  chemical  source  term  has  often  been  considered  to  be  the  main 
problem  in  CFD  simulations  of  turbulent  combustion  (Peters,  2000).  Many 
turbulencechemistry  interaction  models  have  been  developed  and  applied  to  various 
combustion  simulations  since  Spalding  (1971)  proposed  the  eddy-break-up  (EBU)  model 
for  turbulent  premixed  flames.  Most  of  the  models  can  be  placed  in  one  of  the  following 
categories:  eddy-break-up  and  eddy-dissipation  models,  Probability  Density  Function 
(PDF)  models  (including  assumed  PDF  method  and  transport  PDF  method),  laminar 
flamelet  models,  conditional  moment  closures,  and  the  linear  eddy  model  (FEM).  The 
transport  PDF  method  and  the  FEM  model  are  accurate,  but  they  are  also  extremely 
expensive.  Application  of  the  laminar  flamelet  model  is  limited  when  more  than  one  fuel 
stream  is  present.  The  CPU  cost  of  the  assumed  PDF  model  increases  dramatically  when 
applied  to  the  chemical  mechanisms  without  an  analytical  form  of  the  reaction  rate. 

Gran  and  Magnussen  (1996)  introduced  a  new  method  to  implement  finite-rate  chemistry 
into  the  original  eddy  dissipation  concept  (EDC)  model  (Magnussen,  1989),  for  turbulent 
combustion  which  is  an  extension  of  Spalding  s  EBU  model.  Excellent  agreement  with 
experimental  measurements  was  achieved  for  a  bluff-body  stabilized  turbulent  diffusion 
flame.  We  note  that  the  case  geometry  is  similar  to  that  of  a  scramjet  combustor.  In  this 
project,  the  model  was  slightly  modified  and  coupled  with  the  reduced  mechanism.  We 
also  note  that  the  formulation  of  the  model  can  be  directly  coupled  with  ISAT  in  order  to 
reduce  the  CPU  cost. 

The  transport  PDF  method,  where  no  closure  is  needed  for  the  chemical  source  term,  is 
considered  as  the  most  accurate  method  for  turbulent  reactive  flow  simulation.  Due  to  the 
tremendous  amount  of  work  associated  with  implementing  the  PDF  model  into 
VUFCAN  and  limited  time  frame  of  this  project,  the  PDF  method  was  only  applied  using 
a  post-process  approach.  The  density,  velocity,  and  turbulent  kinetic  energy  and 
dissipation  rate  were  interpolated  from  VUFCAN  simulation  results.  The  temperature  and 
species  concentrations  are  recalculated  using  the  PDF  method  without  updating  the 
variables  interpolated  from  VUFCAN. 

Details  of  the  EDC  model  and  PDF  post-processor  are  described  in  the  following 
sections. 

7.1  Eddy  Dissipation  Concept  Model 

The  eddy  dissipation  concept  model  employed  in  this  project  is  an  extension  to  the  EDC 
model  set  forth  by  Gran  and  Magnussen  (1996),  which  is  related  to  the  eddy-break-up 
model  of  Spalding  (1976).  The  rationale  of  the  EDC  model  was  based  on  the  assumption 
that  chemical  reaction  occurs  in  the  regions  where  the  dissipation  of  turbulence  energy 
takes  place.  In  flows  from  moderate  to  intense  turbulence,  these  regions  are  small  and 
isolated,  occupying  only  a  small  fraction  of  the  flow.  These  regions  consist  of  fine 
structures  whose  characteristic  dimensions  are  of  the  order  of  the  Kolmogorov  length 
scale  in  one  or  two  dimensions.  An  empirical  expression  for  mean  reaction  rates  is  then 
given  based  on  this  assumption. 

The  fraction  of  flow  occupied  by  the  isolated  fine  structure  regions  is  expressed  as 
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1/4 


C,F  ,  (6) 

where  k  is  turbulent  kinetic  energy,  is  turbulent  kinetic  energy  dissipation  rate  ,  is  the 
laminar  kinematic  viscosity,  and  CD  2.80  is  a  model  constant  (Gran  and  Magnussen, 

1996).  Note  that  turbulent  kinetic  viscosity  is  t  C  k2  /  .  The  above  formulation  can 
also  be  written  as 

1/4 

CDC  -L  ,  (7) 

t 

where  ,  is  the  laminar  kinetic  viscosity.  This  formula  is  more  convenient  when  the 

turbulence  model  used  in  the  simulation  is  not  the  k-  model.  In  flows  where  turbulent 
viscosity  is  much  larger  than  laminar  viscosity,  the  fraction  of  fine  structures  is  expected 
to  be  small.  Magnussen  (1989)  argued  that  the  volume  fraction  of  the  flow  occupied  by 
the  fine  structures  should  rather  be  modeled  as  2  than  3  because  one  dimension  of  the 
fine  structures  is  comparable  to  large  vortex  structures. 

The  time  scale  for  the  mass  transfer  between  the  fine  structures  and  surroundings  is 
estimated  as 

1/4 

CE  -  ,  (8) 

where  CE  0.167  is  another  model  constant. 

The  mean  reaction  rate  is  then  given  by 

—  (Xi  Y*),  (9) 

where  f 0  is  the  mean  mass  fraction  and  Y*  is  the  mass  fraction  within  the  fine 
structures.  In  this  study,  Y*  is  computed  from 

Y*  o  ( Yt,p,T)dt  ,  (10) 

with  initial  conditions:  Yt  Y .°  at  t  0.  The  pressure  and  enthalpy  are  kept  constant 

during  the  integration.  The  final  form  for  the  mean  reaction  rate  is  written  as 

—  Q  <?iiP,T)dt  If  .  (11) 

From  Equation  (11),  the  turbulent  mean  reaction  rate  given  by  the  EDC  model  is  indeed 
the  averaged  reaction  rate  during  a  time  span,  ,  times  the  volume  fraction  of  the  fine 
structures.  The  relation  between  the  mean  reaction  rate  and  the  parameter,  ,  is  obvious. 
The  life  time  of  the  fine  structures,  ,  is  also  an  important  parameter  in  determining  the 
turbulent  mean  reaction  rate.  The  smaller  the  life  time,  the  closer  the  averaged  reaction 
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rate  calculated  from  the  integration  is  to  the  laminar  reaction  rate.  In  other  words,  smaller 
fine  structure  life  results  in  less  impact  of  turbulence  chemistry  interaction. 

An  Euler  backward  scheme  was  used  to  integrate  Equation  (10),  which  is  given  by 

y;  1  y;  t  (y;  l,p,rnl )  (12) 

where  t  is  the  time  step  of  the  integration,  Y"  and  Y"  1  are  the  species  mass  fractions 

at  time  steps  n  and  n  1 ,  respectively.  The  nonlinear  Equations  (12)  are  solved  using 
Newton’s  method. 

The  integration  procedure  of  Equation  (10)  is  very  computationally  expensive.  The  ISAT 
algorithm,  which  is  discussed  in  Section  6,  was  applied  to  speed  up  the  evaluation  of  the 
mean  reaction  rates.  Note  that  the  tabulation  method  does  not  need  any  change  except  for 
the  inclusion  of  the  parameter,  as  an  independent  variable. 

7.2  PDF  Post-Processor 

Probability  Density  Function  (PDF)  methods  have  been  under  development  for  some 
time  for  modeling  turbulence-chemistry  interactions  with  some  success.  However  CFD 
calculations  with  the  PDF  method  are  too  CPU-intensive  for  modeling  of  practical  3D 
combustor  geometries  with  realistic  chemistry.  In  order  to  assess  the  importance  of 
turbulence-chemistry  interactions,  a  Monte  Carlo  analysis  of  the  joint  scalar  PDF  has 
been  developed  for  use  as  a  post-processing  tool.  The  post-processing  code  reads  in  a 
CFD  solution  containing  the  mean  flow  field  and  turbulence  properties  including 
turbulence  kinetic  energy  and  dissipation  rate.  With  a  given  flow  field,  Monte  Carlo 
simulations  are  performed  to  solve  for  the  joint  scalar  PDF  evolution  equation  using 
fractional  time  steps. 

Figure  18  presents  an  overview  of  the  post-processing  tool.  The  main  routine  PDFMC 
reads  in  the  necessary  information  on  the  turbulent  flows.  As  Monte  Carlo  simulations 
use  notional  particles  to  represent  turbulence  statistics,  numerous  calculations  for  the 
evolution  of  the  chemical  state  are  performed.  The  stiff  kinetics  integrator  is  initialized 
first.  Second,  the  impact  of  turbulent  convection  and  diffusion  is  determined  in 
subroutine  PDFCD  as  the  flow  field  is  assumed  to  remain  unchanged.  Followed  is  the 
main  iteration  routine  STEPMC  which  performs  Monte  Carlo  simulations  for  turbulent 
convection,  diffusion,  molecular  mixing,  and  chemical  reactions  sequentially.  ISAT  is 
used  to  speed  up  kinetics  calculations.  Detailed  input  and  output  files  and  the  associated 
manipulations  by  the  post-processing  tool  are  presented  in  Figure  19. 
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Figure  18.  Schematic  of  PDF  Post-Processing  Tool  Showing  that  Results  from  CFD 
are  Read  in  and  the  Main  Routine  Starts  the  Monte  Carlo  Simulations. 


Input  files: 

Flow  field,  boundary  conditions 
geometry,  fuel,  oxidizer  properties 


PDF  -post  processing 


Output  files: 

PDF  restart  file,  contour  plots 


Plotting  processing  codes 


PDF  scatter  plots 
Tecplot  contour  plots 


Figure  19.  Schematic  of  Details  of  PDF  Post-Processing  Tool  Showing  Various 
Input  and  Output  Files  Associated  with  Different  Parts  of  the  Process. 
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8.0  Implementation  of  New  Techniques  in  VULCAN 

8.1  VULCAN  Overview 

VULCAN  is  a  Navier-Stokes  flow  solver  that  is  currently  maintained  and  distributed 
(restricted  to  the  U.S.)  by  the  Hypersonic  Airbreathing  Propulsion  Branch  of  NASA 
Langley  Research  Center.  For  details,  see  the  Vulcan  home  page  (http://vulcan- 
cfd.larc.nasa.gov/)  (White  &  Morrison,  1999). 

VULCAN  is  a  multigrid,  multiblock,  structured,  finite-volume  code,  developed  for 
solving  the  spatially  elliptic  and  parabolized  forms  of  the  equations  governing  3-D, 
turbulent,  calorically  perfect  and  nonequilibrium  chemically  reacting  flows.  Space 
marching  algorithms  developed  to  improve  convergence  and/or  reduce  computational 
cost  were  implemented  as  well  as  elliptic  methods  for  solving  flows  with  large  regions  of 
subsonic  flow.  A  full  approximate  storage,  full  multigrid  scheme  was  also  implemented 
to  accelerate  convergence  of  either  the  elliptic  or  space  marching  schemes. 
Compressibility  corrected  forms  of  the  k-  and  k-  ,  two  equation  turbulence  models, 
were  implemented  that  are  suitable  for  high-speed  flows.  In  addition,  a  compressible, 
pressure  gradient  corrected  turbulent  law-of-the-wall  matching  function  was  implemented 
that  decreased  wall  grid  spacing  sensitivity.  Turbulence-chemistry  interaction  models 
were  also  implemented  using  the  assumed  PDF  methods. 

For  the  simulations  completed  in  this  work,  turbulence  was  modeled  with  the  Menter 
(1992)  baseline  (BSL)  two-equation  turbulence  model.  The  Menter  BSL  model  couples 
the  standard  forms  of  the  k-  and  k-  models.  The  BSL  model  invokes  the  k-  model 
near  solid  surfaces  and  smoothly  transitions  to  the  k-  model  in  the  outer  portion  of  the 
boundary  layer  and  in  regions  of  free  shear.  At  solid  surfaces,  the  wall  matching 
procedure  of  Wilcox  (1989)  was  employed.  Also,  the  compressibility  correction  of 
Wilcox  (1993)  was  used  to  model  the  reduction  in  mixing  associated  with  high 
convective  Mach  numbers. 


8.2  Implementation  of  Reduced  Mechanisms 

To  implement  reduced  mechanisms  in  the  VULCAN  code,  a  numerical  Jacobian 
capability  was  first  implemented  in  VULCAN.  This  was  necessary  because  the  implicit 
treatment  of  chemical  source  terms  in  VULCAN  requires  the  Jacobian  matrices  of  the 
reaction  rates,  allowing  VULCAN  to  use  the  (nonanalytical)  chemical  rate  source  terms 
from  CARM-formulated  reduced  mechanisms.  CARM  produces  a  subroutine  that 
calculates  chemical  source  terms  for  a  given  set  of  conditions,  rather  than  terms  which  fit 
an  analytical  form  such  as  Arrhenius  rate  coefficients. 

The  numerical  Jacobian  has  been  implemented  into  VULCAN  including  both  the  species 
derivative  and  temperature  derivative  contributions.  It  is  calculated  using  the  formula 


/  ) 


(13) 
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where  is  given  by: 


10  7  10  6 


The  implementation  has  been  checked  with  two  hand  calculations  using  a  conventional 
seven-specie,  seven-reaction  Arrhenius  hydrogen-air  model.  The  numerical  values  agreed 
to  10  significant  figures  for  the  species  derivative  and  to  3  significant  figures  for  the 
temperature  derivative.  Also,  the  terms  of  the  numerical  Jacobian  have  been  compared  to 
the  corresponding  terms  of  the  analytical  Jacobian  as  implemented  into  VULCAN.  Good 
agreement  was  found  between  the  terms  of  the  analytical  Jacobian  in  VULCAN  and  the 
new  numerical  Jacobian. 

The  implementation  was  then  tested  for  an  ignition  delay  case  (  =  0.25,  T  =  1000  K,  P  = 
1  atm)  using  the  seven-specie  seven-reaction  Arrhenius-form  hydrogen-air  model. 
Ignition  delay  was  taken  as  the  time  for  the  temperature  to  increase  400  K.  The  ignition 
delay  calculated  using  Chemkin  and  the  seven-by-seven  hydrogen-air  model  was  1.580E- 
04  s.  In  VULCAN,  a  three-stage  Runge-Kutta  method  was  used.  The  ignition  delay 
values  obtained  with  analytical  and  numerical  Jacobians  were  1.566e-4  s  and  1.564e-4  s, 
respectively.  Figure  20  shows  the  time  history  of  temperature  for  the  calculations.  The 
solutions  with  the  two  methods  are  very  close. 

With  the  numerical  Jacobian  capability  added  to  VULCAN,  the  implementation  of  the 
reduced  mechanism  is  straightforward.  The  gas  temperature,  pressure  and  species  mass 
fractions  are  passed  to  the  reduced  mechanism  subroutine,  CKWYP,  and  the  chemical 
source  terms  are  returned. 
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Figure  20.  Temperature  Versus  Time  Using  Analytical  Jacobian  and  Numerical 
Jacobian  for  an  Ignition  Delay  Calculation. 

8.3  Implementation  of  ISAT 

Figure  21  shows  the  overview  of  the  interconnection  between  a  reactive  flow  code,  the 
ISAT  algorithm,  and  the  reduced  mechanism.  The  CFD  code  passes  the  species  mass 
fractions,  temperature,  and  pressure  to  the  ISAT  routine.  The  ISAT  routine  then  traverses 
the  existing  binary  tree  table.  If  an  entry  within  the  ellipsoidal  of  accuracy  is  found,  an 
interpolation  is  performed  and  the  interpolated  reaction  rates  and  the  Jacobian  matrix 
used  for  interpolation  are  returned  to  VULCAN.  If  such  an  entry  is  not  found,  the  ISAT 
routine  then  calls  CKWYP  to  calculate  the  reaction  rates  and  Jacobian.  A  new  entry  is 
added  to  the  table,  and  the  reaction  rates  and  Jacobian  are  returned  to  VULCAN.  The 
implicit  method  used  in  VULCAN  requires  the  evaluation  of  the  Jacobian  matrix  of  the 
chemical  source  term,  which  is  available  in  the  ISAT  table.  It  thus  saves  significant  CPU 
time  for  the  ISAT  subroutine  to  return  not  only  the  reaction  rate  vector,  but  the  Jacobian 
matrix  required  by  VULCAN.  The  only  additional  cost  introduced  by  ISAT  is  traversing 
the  binary  tree  data  structure  to  search  for  a  close  entry  to  the  inquiry  vector.  The  cost  of 
this  operation  is  very  small  compared  to  direct  evaluation  of  the  reaction  rates. 

The  memory  size,  M,  occupied  by  the  ISAT  table  can  be  estimated  as  M=  8 N(n  +6/7+ 1 4) 
bytes,  where  n  is  the  number  of  species,  and  N  is  the  number  of  table  entries.  For  a  20- 
specie  mechanism  and  a  table  size  of  50,000,  the  memory  size  occupied  by  the  table  is 
approximately  214  MB.  It  was  found  that  the  ISAT  table  size  increases  very  fast  in  the 
early  stage  of  the  simulations  due  to  the  fast  change  in  species  concentrations  and 
temperature.  Many  table  entries  created  early  are  rarely  used  in  the  later  iterations  or  time 
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steps.  A  tree-trimming  method  was  developed  to  control  the  table  size.  An  upper  limit  is 
set  for  the  size  of  binary  tree.  Once  the  tree  is  full,  the  table  entries  that  have  never  been 
used  during  the  most  recent  few  time  steps  are  trimmed  from  the  table. 
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Figure  21.  Overview  of  the  Interconnection  between  a  Reactive  Flow  Code,  ISAT, 
and  a  Reduced  Mechanism. 

8.4  Implementation  of  EDC  Model 

Implementation  of  the  EDC  model  in  VULCAN  is  very  similar  to  that  of  the  reduced 
mechanism,  as  shown  in  Figure  22.  The  target  functions  of  ISAT  are  the  turbulent  mean 
reaction  rates,  and  the  independent  variables  for  ISAT  are  species  mass  fractions, 
pressure,  temperature,  and  the  life  time  of  the  fine  structures.  Utilizing  the  EDC  model 
without  ISAT  is  very  CPU  intensive,  due  to  the  large  number  of  calls  to  the  reduced 
mechanism  routine  in  order  to  integrate  the  mean  reaction  rates.  Thus,  the  efficiency  of 
ISAT  is  expected  to  be  significantly  better  than  when  using  laminar  chemistry. 


Figure  22.  Overview  of  the  interconnection  between  a  reactive  flow  code,  ISAT, 
EDC  model,  and  a  reduced  mechanism. 
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9.0  Evaluation  of  New  Techniques  in  VULCAN 

Initial  evaluations  of  the  reduced  mechanism  and  ISAT  implementations  into  VULCAN 
were  performed  for  2-D  supersonic  reacting  flows.  This  allowed  debugging  and 
examination  of  stability,  accuracy,  and  parallelization  issues  to  be  performed  on  cases 
with  comparatively  small  CPU  demands.  These  reactions  are  discussed  in  Section  9. 1 . 
Simulations  were  completed  to  investigate  impacts  of  turbulence-chemistry  interactions 
on  calculation  of  the  mean  chemical  reaction  rates.  A  discussion  of  these  results  is 
provided  in  Section  9.2. 

9.1  Evaluation  of  Reduced  Mechanisms  with  ISAT 

The  case  considered  is  the  coaxial  jet  flow  experiment  conducted  by  Evans  et  al.  (1978). 
The  experiment  consisted  of  coaxial  injection  of  a  cold,  Mach  2,  hydrogen  jet  at  matched 
pressure  into  a  hot,  Mach  1.9,  vitiated  air  stream.  Figure  23  is  a  schematic  of  this 
experiment. 


Figure  23.  Schematic  of  the  Supersonic  Diffusion  Flame  Experiment  of  Evans  et  al. 
(1978). 

The  fuel  nozzle  was  conical  with  a  5°  exit  half-angle.  The  computational  domain 
extended  seven  jet  diameters  ( d  =  0.009525  m)  upstream  of  the  nozzle  exit  to  capture  the 
flow  angularity  at  the  inner  nozzle  exit  and  to  include  a  boundary  layer  in  the  exit 
profdes.  Modeling  the  flow  upstream  of  the  jet  exit  was  found  to  yield  better  agreement 
with  the  experimental  data.  The  computational  domain  also  extended  30  diameters 
downstream  of  the  nozzle  exit  and  2  diameters  in  the  transverse  direction.  The  modeled 
flow  conditions  are  given  in  Table  3.  The  computational  grid  downstream  of  the  nozzle 
exit  plane  that  was  used  for  the  majority  of  the  calculations  was  discretized  with  120  and 
68  cells  in  the  streamwise  and  transverse  directions,  respectively.  Calculations  were 
performed  with  grids  containing  240x136  and  60x34  cells  to  assess  grid  sensitivity  in  the 
results.  Complete  grid  independence  was  not  obtained,  but  differences  between  the 
solutions  on  the  medium  and  fine  grids  were  minor.  The  results  shown  in  the  succeeding 
figures  were  obtained  on  the  medium  grid  (120x68). 
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Table  3.  Conditions  for  Supersonic  Jet  Flame  Simulations. 


Specified  condition 

Fuel  jet 

Outer  jet 

Mach  number 

— 

1.9 

Temperature  (K) 

— 

1495 

Pressure,  p  (bar) 

— 

1.0 

Total  temperature  (K) 

452 

— 

Total  pressure  (bar) 

8.6068 

— 

02mass  frac. 

0.0 

0.241 

N2mass  frac. 

0.0 

0.478 

H20  mass  frac. 

0.0 

0.281 

Fuel  mass  frac. 

1.0 

0.0 

Solutions  were  obtained  in  three  steps.  First,  flow  upstream  of  the  jet  nozzle  exit  plane 
was  calculated.  The  resulting  solution  was  used  as  the  inflow  boundary  condition  for  the 
flow  calculation  downstream  of  nozzle  exit  plane.  A  no-slip,  adiabatic  boundary 
condition  was  applied  at  the  injector  lip,  a  symmetry  boundary  condition  was  applied 
along  the  centerline,  and  zeroth  order  extrapolation  was  used  at  the  outflow  boundary  and 
at  the  upper  domain  boundary.  Secondly,  the  coaxial  jet  flow  field  calculation  was 
advanced  with  reactions  disabled.  Finally,  the  calculation  was  continued  with  the 
reactions  enabled.  This  same  procedure  was  followed  for  both  the  CARM-produced  and 
global  Arrhenius-form  models.  Similar  convergence  behavior  was  observed  in  both 
cases. 

9.1.1  Supersonic  Hydrogen  Flame  Simulations 

Solutions  were  obtained  for  both  reduced  and  detailed  H2-air  kinetics  models.  The 
detailed  mechanism  was  a  subset  of  the  GRI  3.0  mechanism  (GRI  mech  web  page)  and 
included  the  nine  species  H2,  O2,  H2O,  OH,  O,  H,  HO2,  H2O2,  and  N2.  All  of  the  reactions 
involving  these  species  from  GRI  3.0  were  included  except  one,  which  was  omitted 
because  pressure-dependent  reactions  are  currently  not  supported  by  VULCAN.  The 
resulting  mechanism  had  27  reactions. 

For  comparison  purposes,  two  mechanisms  were  evaluated.  The  first  is  the  CARM- 
produced  reduced  kinetic  mechanism  using  the  seven  species  H2,  O2,  H2O,  O,  H,  HO2, 
and  N2,  which  was  described  earlier.  The  second  is  a  global  mechanism  obtained  by 
eliminating  HO2  and  H2O2  from  the  mechanism  of  Jachimowski  (1988),  which  yields  a 
seven-step,  seven-species  Arrhenius-form  mechanism  which  we  refer  to  in  the  following 
as  the  seven-by-seven  mechanism. 

Contours  of  the  Mach  number,  temperature,  fuel-air  E.R.,  and  the  ratio  of  turbulent-to- 
laminar  viscosity  from  the  solution  with  the  CARM  seven-species  reduced  mechanism 
are  shown  in  Figure  24.  The  hydrogen  jet  expanded  from  sonic  conditions  to  a  nominal 
Mach  number  of  2.0  at  the  nozzle  exit.  Further  expansion  elevated  the  Mach  number  to  a 
maximum  of  2.72.  The  hydrogen  jet  turned  back  toward  the  centerline  at  x/D  ~  1.0, 
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recompressing  the  jet.  The  minimum  Mach  number  at  the  exit  of  the  domain  (x/D  =  30) 
was  1.26. 

A  mixing  layer  formed  in  the  wake  behind  the  nozzle  lip  between  the  hydrogen  jet  and 
the  vitiated  air  stream.  The  wake  also  generated  relatively  high  levels  of  turbulent  kinetic 
energy,  and  thus,  high  eddy  viscosity  values,  which  were  convected  downstream.  The 
fuel-air  E.R.s  indicate  a  thin  surface  of  approximately  stoichiometric  proportions  along 
the  outer  edge  of  the  shear  layer  that  formed  between  the  two  streams.  Ignition  was 
observed  at  x/D  ~  2.  The  flame  occurred  near  the  stoichiometric  surface.  Contours  of  the 
02,  N2,  H20,  and  H2  mass  fractions  from  the  same  solution  are  shown  in  Figure  25. 
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Figure  24.  Contours  of  Mach  Number,  Fuel-Air  E.R.,  Ratio  of  Turbulent  to 
Laminar  Viscosity  and  Temperature  for  the  Solution  Obtained  with  the 
CARM  Seven-Species  Reduced  Chemical  Kinetic  Mechanism. 
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Figure  25.  Contours  of  Mass  Fraction  of  O2,  N2,  H20,  and  H2  for  the  Solution 
Obtained  with  the  CARM  Seven-Species  Reduced  Chemical  Kinetic  Mechanism. 

Oxygen  penetrated  into  the  fuel-rich  region  upstream  of  the  flame  surface.  This  oxygen 
diffused  into  the  fuel  jet  and  reacted  with  the  fuel.  Downstream  of  ignition,  little  oxygen 
penetrated  through  the  flame  zone.  N2,  H20,  and  H2  diffuse  in  the  region  below  the  flame 
surface.  Contours  of  H20  reveal  that  little  diffusion  occurs  above  the  shear  layer  due  to 
low  values  of  turbulent  viscosity  at  the  outer  edge  of  the  shear  layer. 

Profiles  of  the  mass  fraction  of  the  major  species  02,  N2,  H20,  and  H2  are  compared  for 
the  solutions  with  three  kinetics  models  and  the  experimental  data  in  Figure  26. 
Experimental  data  for  other  species  and  temperature  were  not  available.  As  reflected  in 
the  N2  profiles,  near  the  centerline,  the  shear  layer  mixing  is  slightly  underpredicted  at  the 
first  measurement  station  but  predicted  well  at  the  succeeding  three  stations.  However, 
the  extent  of  mixing  near  the  edge  of  the  shear  layer  (x/d  ~  0.7)  is  greater  in  the 
calculations  than  in  the  experimental  data.  The  measured  values  for  02  and  H20  in  the 
vitiated  stream  (x/d  >  1.2)  are  observed  to  deviate  from  their  nominal  values.  As  seen  in 
Figure  25,  oxygen  penetrated  into  the  fuel  rich  region,  and  positive  values  of  oxygen 
mass  fraction  are  predicted  below  the  flame  at  x/d  =  8.26.  Each  of  the  kinetics  models 
predicted  more  oxygen  than  was  observed  in  the  measurements,  with  the  CARM  reduced 
model  predicting  the  least  oxygen  and  the  seven-by-seven  model  the  most.  Both  the 
CARM  reduced  mechanism  and  the  detailed  model  predicted  negligible  02  at  the  last 
three  stations  in  accord  with  the  measurements,  while  the  seven-by-seven  model 
predicted  significant  levels  of  02.  This  demonstrates  the  presence  of  kinetic  effects 
retained  in  the  CARM  model  and  detailed  chemistry,  but  absent  from  the  seven-by-seven 
mechanism. 
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Figure  26.  Comparison  of  Species  Profiles  between  Calculation  and  Experiment. 
9.1.2  Supersonic  Ethylene  Flame  Simulations 

Calculations  using  the  same  geometry  and  boundary  conditions  were  repeated  with 
ethylene  (C2H4)  as  the  fuel  using  the  20-species  reduced  mechanism  described 
previously.  The  total  pressure  and  temperature  of  the  C2H4  jet  were  the  same  as  for  the  H2 
jet  to  again  produce  matched  pressure  conditions.  No  numerical  stability  problems  were 
encountered.  Shown  in  Figure  27  are  contours  of  the  temperature  and  the  mass  fractions 
of  O2,  N2,  CO2,  and  C2H4.  The  standoff  flame  distance  has  increased  from  x/d  ~  2.0  to  x/d 
~  in  the  hydrogen-air  flame  and  x/d  ~  17.0  in  the  ethylene-air  flame.  The  calculation 
demonstrates  the  capability  for  modeling  high-speed  combustion  problems  using 
hydrocarbon  fuels  with  CARM-produced  kinetic  mechanisms  in  a  state-of-the-art  CFD 
solver.  Unfortunately,  very  little  experimental  data  is  available  for  supersonic 
hydrocarbon  flames. 
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Figure  27.  Contours  of  Temperature  and  Mass  Fraction  of  O2,  N2,  CO2,  and  C2H4 
for  the  Solution  Obtained  with  the  20-Species  C2H4  Reduced  Chemical  Kinetic 
Mechanism. 

The  error  tolerance  for  ISAT  for  both  cases  was  set  to  0.02,  with  a  maximum  table  size  of 
50,000  entries.  After  several  hundred  iterations,  the  table  size  stabilized  at  around  30,000 
entries.  The  number  of  new  entries  added  to  the  table  at  each  time  step  was  less  than  10. 
The  average  CPU  time  per  time  step  for  the  mechanism  is  listed  in  Table  4.  For  the 
hydrogen  flame,  the  CPU  cost  of  the  reacting  flow  simulation  with  directly  evaluated 
reaction  rates  was  only  1.8  times  that  of  the  nonreacting  case.  Still,  a  speedup  factor  of 
1.34  was  achieved  with  ISAT.  For  the  ethylene  flame  case,  the  CPU  cost  of  the  reacting 
flow  simulation  without  ISAT  is  about  82  times  that  of  nonreacting  flow.  Numerically 
evaluating  the  reaction  rate  Jacobian  matrix  is  extremely  expensive  for  the  CARM- 
produced  reduced  mechanisms.  More  than  98%  of  the  total  CPU  time  was  consumed  by 
the  reduced  mechanism  subroutine.  Using  ISAT,  the  total  CPU  time  was  cut  by  a  factor 
of  36,  making  the  reacting  flow  comparable  in  CPU  time  to  the  nonreacting  flow 
simulation. 


Table  4.  Average  CPU  Time  for  a  Single  Time  Step  for  the  Jet  Flame  Simulations. 


Mechanism 

7-specie  H2 

20-specie  C2H4 

CPU,  nonreacting 

0.35  s 

0.93  s 

CPU,  reacting,  direct 

0.63  s 

76.4  s 

CPU,  reacting,  ISAT 

0.47  s 

1.8  s 

ISAT  speedup  factor 

1.34 

36.4 
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Calculations  of  the  Evans  et  al.  (1978)  diffusion  flame  geometry  case  have  been  repeated 
for  comparison  to  other  ethylene  kinetic  models.  A  solution  was  obtained  with  a  six 
species,  three-reaction  ethylene  mechanism  which  has  been  used  previously  in  scramjet 
combustor  simulations  (Baurle  &  Eklund,  2002).  The  CPU  time  compared  to  the  20- 
species  ethylene  reduced  kinetics  model  with  and  without  the  use  of  ISAT  is  shown  in 
Table  5.  The  timings  were  obtained  on  a  single  Compaq  SC-45  processor  for  each  of  the 
three  calculations.  ISAT  reduced  the  CPU  time  by  a  factor  of  42. 

Table  5.  CPU  Times  for  Calculations  Performed  on  One  Processor. 


Mechanism 

Seconds  per  iteration 

C2H4  (6  by  3) 

0.54 

C2H4  (reduced  20  specie) 

76.4 

C2H4  (reduced  20  specie)  +ISAT 

1.8 

9.1.3  Parallelization 

Further  calculations  were  performed  using  the  20-species  C2H4  reduced  mechanism  and 
ISAT,  but  using  multiple  processors  on  a  Compaq  SC-45.  The  results  of  the  timings  are 
shown  in  Table  6.  The  relatively  poor  speedup  is  due  to  the  small  grid  employed  (121  by 
69)  and  the  consequent  large  relative  ratio  of  communication  to  computation,  which 
increases  as  the  number  of  processors  increases. 

Table  6.  CPU  Times  for  2-D  Ethylene  Jet  Flame  Calculations  (with  ISAT)  on 
Multiple  Processors. _ _ _ 


Processors 

Seconds  per  iteration 

Speedup 

1 

1.8 

— 

2 

1.1 

1.6 

4 

0.7 

2.6 

8 

0.5 

3.5 

12 

0.5 

3.6 

24 

0.8 

2.3 

No  convergence  problems  were  encountered  using  multiple  processors,  as  shown  in 
Figure  28,  which  shows  the  residual  histories  for  the  solutions  obtained  using  the  reduced 
kinetics  model  with  ISAT  and  various  number  of  processors.  The  solution  labeled 
DIRECT  refers  to  the  solution  obtained  without  ISAT. 
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Figure  28.  Residual  History  for  Solutions  obtained  with  ISAT  and  without  ISAT 
(DIRECT). 

Discrepancies  were  initially  observed  among  solutions  obtained  using  different  numbers 
of  processors.  Figure  29  displays  contours  of  OH  mass  fraction  for  the  solution  obtained 
with  the  reduced  kinetics  model  without  ISAT,  while  Figures  30  and  31  show  contours  of 
OH  mass  fraction  in  the  vicinity  of  the  flame  region  for  the  solutions  using  ISAT.  The 
solution  with  ISAT  and  one  processor  (red  contours)  is  indistinguishable  from  the 
solution  without  ISAT  (black  contours);  however,  the  solutions  with  multiple  processors 
are  clearly  different  from  each  other. 
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Figure  29.  Contours  of  OH  Mass  Fraction  Obtained  with  the  Reduced  Kinetics 
Model  without  ISAT. 
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Figure  30.  Contours  of  OH  Mass  Fraction  Obtained  with  ISAT  Using  1  and  2 
Processors  and  without  ISAT  (DIRECT).  See  Figure  28  for  the  color  key. 
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Figure  31.  Contours  of  OH  Mass  Fraction  Obtained  with  ISAT  using  4,  8  12  and  24 
Processors.  See  Figure  28  for  the  color  key. 

Differences  in  the  solutions  are  more  clearly  apparent  in  Figures  32  and  33,  which  show 
profiles  of  OH  mass  fraction  at  two  axial  locations.  The  values  of  OH  mass  fraction 
above  the  flame  region  (y  >  0.013  m  in  Figure  30  andy  >  0.017  m  in  Figure  31)  are 
nonzero  for  the  solutions  obtained  with  multiple  processors.  The  calculation  was  repeated 
with  four  processors  and  with  ISAT  the  error  tolerance  tightened  from  0.05  to  0.01. 
Shown  in  Figure  34  is  the  residual  history.  Figure  35  displays  the  contours  of  OH  mass 
fraction  and  Figures  36and  37  show  the  OH  mass  fraction  profiles.  The  solution  with  the 
tighter  ISAT  error  tolerance  is  now  indistinguishable  from  the  solution  obtained  without 
ISAT. 
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Figure  32.  OH  Mass  Fraction  Profiles  at  x  =  0.182  m. 


Figure  33.  OH  Mass  Fraction  Profiles  at  x  =  0.286  m  (exit  of  the  domain). 


48 


Figure  34.  Residual  History  with  the  ISAT  Error  Tolerance  Reduced  to  0.01. 


Figure  35.  Contours  of  OH  Mass  Fraction  obtained  with  ISAT  and  the  Error 
Tolerance  Set  to  0.05  and  0.01  and  Without  ISAT  (DIRECT).  See  Figure  34  for  the 
color  key. 
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Figure  36.  OH  Mass  Fraction  Profiles  at  x  =  0.182  m  with  the  Error  Tolerance 
Reduced  to  0.  01. 


Figure  37.  OH  Mass  Fraction  Profiles  at  x  =  0.286  m  with  the  Error  Tolerance 
Reduced  to  0.  01. 
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It  is  believed  that  the  solution  differences  seen  using  ISAT  and  multiple  processors  are 
due  to  separate  ISAT  tables  that  are  generated  for  each  processor.  Differences  among  the 
tables  will  require  a  tighter  overall  tolerance  on  the  ISAT  interpolation  to  achieve  the 
same  results.  Reducing  the  error  tolerance  caused  the  CPU  time  to  increase  from  0.7 
seconds/iteration  to  1 .2  seconds/iteration. 

Finally,  a  3-D  grid  was  used  to  further  evaluate  the  parallel  speedup  of  the  ISAT  routine. 
The  grid  was  obtained  from  the  axisymmetric  grid  by  sweeping  180°  in  36  intervals;  the 
resulting  grid  size  was  121  by  69  by  37  cells.  Shown  in  Figure  38  is  a  cross  section  of  the 
3-D  grid.  The  CPU  timings  for  the  ISAT  routine  on  the  3-D  grid  are  given  in  Table  7. 
Note  that  for  the  parallel  computations  both  the  axisymmetric  and  3-D  grids  were  split  in 
the  axial  direction. 


Y 

Figure  38.  Outflow  Plane  for  the  3-D  Grid. 


Table  7.  CPU  Times  for  Calculations  Performed  on  Multiple  Processors  Using  a  3- 
D  Grid. 


Processors 

Sec.  per  iteration 

Speedup 

Efficiency 

1 

73.6 

— 

1.0 

2 

41.1 

1.8 

0.90 

4 

25.1 

2.9 

0.73 

8 

16.2 

4.5 

0.57 

12 

11.3 

6.5 

0.54 

24 

8.63 

8.5 

0.36 
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9.1.4  2-D  A-Heptane  Case 

The  19-species  H-heptane  reduced  mechanism  has  been  successfully  applied  to  simulate 
the  axisymmetric  diffusion  flame  using  VULCAN.  The  inflow  conditions  of  the  flame 
simulation  are  similar  to  the  experimental  setup  of  Evans  et  al.  (1978)  except  that  the 
static  temperature  of  the  air  jet  was  again  increased  by  200  K  to  facilitate  the  ignition  of 
the  flame.  The  grid  size  for  the  simulation  is  144  by  69  cells.  Shown  in  Figure  39  are  the 
contour  plots  of  static  temperature  and  mass  fractions  of  the  species  C7H16,  O2,  and  H. 
The  speedup  factor  when  using  ISAT  for  this  case  is  approximately  10. 


Figure  39.  Contour  Plots  of  Static  Temperature  and  C7H16,  O2,  and  H  Mass 
Fractions  for  the  2-D  Heptane  Diffusion  Flame. 

9.1.5  2-D  JP-8  Case 

The  3  5 -species  JP-8  reduced  mechanism  has  been  successfully  applied  to  simulate  the 
axisymmetric  diffusion  flame  with  VULCAN.  The  inflow  conditions  of  the  flame 
simulation  are  again  similar  to  the  experimental  setup  of  Evans  et  al.  (1978)  except  that 
the  static  temperature  of  the  air  jet  was  increased  by  200  K  to  facilitate  the  ignition  of  the 
flame.  The  grid  size  for  the  simulation  is  144  by  69  cells.  Shown  in  Figure  40  are  contour 
plots  of  static  temperature,  Mach  number,  and  mass  fractions  of  the  six  species  in  the  fuel 
surrogate,  O2,  CO,  and  CO2.  This  case  demonstrates  the  capability,  developed  during  this 
project,  of  taking  a  large,  complex  mechanism  and  fuel  surrogate  from  the  literature  and 
creating  a  reduced  mechanism  that  can  be  incorporated  into  CFD  simulations. 
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Figure  40.  Contour  Plots  of  Static  Temperature,  Mach  Number,  and  Species  Mass 
Fractions  for  the  2-D  JP-8  Diffusion  Flame. 

9.1.6  ISAT  Impacts  for  Complex  Fuels 

The  CPU  times  of  simulations  of  the  /7-heptane  and  JP-8  flames  are  summarized  in  Table 
8.  Note  that  the  computer  platform  used  to  run  the  two  cases  is  different  from  that  for  the 
H2  and  C2H4  case  shown  earlier  in  this  section.  For  comparison  purposes,  the  CPU  times 
for  C2H4  case  are  also  listed  in  the  table.  The  error  tolerance  for  each  of  the  three  cases 
was  2%.  The  speedup  factors  of  ISAT  for  the  77-heptane  and  JP-8  flames  are  1 1  and  12, 
respectively.  Both  are  much  smaller  than  the  factor  of  36  achieved  in  the  C2H4  case.  The 
ratio  of  CPU  times  between  the  reacting  case  and  the  mixing  case,  which  can  be  used  as 
an  indicator  of  how  expensive  the  mechanism  is,  is  72,  15,  and  34,  respectively,  for  the 
three  mechanisms.  A  smaller  speedup  factor  is  expected  for  the  relatively  inexpensive  n- 
heptane  mechanism.  In  VULCAN,  the  implicit  method  requires  the  inversion  of  the 
Jacobian  matrix,  the  chemical  source  term,  of  which  the  computational  cost  is 
proportional  to  77 3 ,  where  n  is  number  of  species.  In  the  35-specie  JP-8  mechanism,  the 
CPU  cost  for  the  inversion  is  about  5  times  that  of  the  19-specie  /7-heptane  mechanism. 
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This  cost  remains  the  same  with  and  without  ISAT.  It  results  in  a  lower  speedup  for  the 
JP-8  flame  case.  The  exact  measurement  of  the  CPU  time  consumed  by  different 
procedures  in  VULCAN,  however,  is  not  available  at  this  time.  The  CPU  time  of  the 
reacting  case  for  each  of  the  mechanisms,  is  within  3  times  the  CPU  time  of  the  mixing 
case. 


Table  8.  Average  CPU  Times  for  a  Single  Time  Step  of  Simulations  of  C2H4,  n- 
Heptane,  and  JP-8  Flames. 


Mechanisms 

M  ixing 

Reacting 

w/o  ISAT 

Reacting 

w/  ISAT 

Speedup 

C2H4,  20-specie 

1.43 

102.9 

2.82 

36 

N -Heptane,  19- species 

1.32 

19.6 

1.72 

11 

JP-8,  3 5- species 

2.11 

72.1 

5.94 

12 

9.1.7  Summary  of  2-D  Implementation 

The  2-D  cases  described  in  this  section  demonstrate  that  detailed  chemical  kinetic 
mechanisms  can  be  reduced  using  CARM  and  successfully  integrated  into  VULCAN. 
Use  of  these  reduced  mechanisms  has  been  demonstrated  to  include  kinetic  information 
that  is  lost  when  using  comparable  simplified  Arrhenius-form  models.  Reduced 
mechanisms  for  larger  hydrocarbons  can  be  very  slow  to  evaluate,  but  the  use  of  ISAT 
makes  the  CPU  time  needed  for  CARM  reduced  mechanisms  comparable  to  that  of 
similarly  sized  global  kinetics.  The  speedup  for  multiple  processors  when  using  ISAT  is 
not  as  large  as  anticipated,  but  the  speedup  improves  for  3-D  cases.  2-D  cases  that  were 
run  using  reduced  mechanisms  for  /7-heptane  and  JP-8  demonstrate  the  feasibility  of 
creating  and  implementing  reduced  kinetics  for  large  hydrocarbons  into  VULCAN. 
Unfortunately,  no  experimental  data  are  available  for  comparison. 

9.2  Turbulence-Chemistry  Interaction 
9.2.1  EDC  Model 

The  EDC  turbulent  combustion  model  was  applied  to  simulate  the  hydrogen  and  ethylene 
flames  described  in  Sections  9.1.1  and  9.1.2.  The  grid  size  for  all  the  cases  shown  here  is 
61  by  35  cells.  For  the  hydrogen  flame  simulations,  Figure  41  compares  the  temperature 
and  H2  mass  fraction  fields  predicted  with  the  EDC  model  and  laminar  chemistry.  The 
flame  structures  look  very  similar  for  both  cases.  The  flame  is  detached  in  both  cases. 

The  liftoff  height  of  the  flame  predicted  with  the  EDC  model  is  slightly  larger  than  that 
with  laminar  chemistry.  This  implies  that  the  turbulence-chemistry  interaction  might  be 
important  when  the  flame  is  close  to  blowout  conditions. 
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Figure  41.  Temperature  and  H2  Mass  Fraction  Fields  Predicted  with  EDC  Model 
and  Laminar  Chemistry. 


The  H2  mass  fraction  along  the  centerline  is  shown  in  Figure  42.  The  potential  core  of  the 
EDC  flame  is  slightly  shorter  than  that  of  laminar  chemistry,  and  the  downstream 
centerline  H2  mass  fraction  is  slightly  lower  than  that  of  laminar  chemistry  downstream. 
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Figure  42.  Distribution  of  H2  Mass  Fraction  Along  the  Centerline. 

Figures  43  and  44  show  the  radial  profiles  of  FL  and  O2  mass  fractions  at  two  axial 
locations,  respectively.  Again  one  sees  very  similar  flame  structure.  The  larger  liftoff 
height  of  the  EDC  model  allows  more  air  mixed  into  the  fuel  stream,  which  results  in 
higher  O2  concentration  near  the  centerline. 
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Figure  43.  Radial  Distribution  of  H2  and  O2  Mass  Fractions  at  x=8.26D. 
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Figure  44.  Radial  Distribution  of  H2  and  O2  Mass  Fractions  at  x=15.5D. 

In  Figure  45,  two  important  parameters  in  the  EDC  model,  the  life  time  of  fine  structures, 
,  and  the  fraction  of  the  regions  occupied  by  the  fine  structures,  are  shown.  The 
parameter  is  smaller  inside  the  core  of  the  jet,  and  larger  outside.  Throughout  the 
computational  domain,  it  is  on  the  order  of  10  5  seconds,  which  is  much  smaller  than  the 
typical  ignition  delay  time  of  the  hydrogen.  As  discussed  in  Section  7.1,  the  averaged 
reaction  rate  from  Equation  (1 1)  is  expected  to  be  close  to  the  laminar  reaction  rate.  The 
small  values  of  the  life  time  of  the  fine  structures,  which  resulted  from  the  high-speed 
nature  of  the  jet,  are  insignificant  in  their  contribution  to  the  mean  reaction  rate.  The 
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fraction  of  the  fine  structure,  ,  is  close  to  1  outside  the  jet  shear  layer,  where  turbulence 
is  not  very  strong.  Just  downstream  of  the  hydrogen  duct,  there  is  a  region  where  is  very 
small.  This  results  in  a  small  turbulent  mean  reaction  rate,  and  delayed  ignition  of  the 
flame. 
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Figure  45.  Contour  Plots  of  EDC  Model  Parameters,  (top,  in  seconds)  and 
(bottom). 


Comparisons  of  predicted  temperature  and  C2H4  mass  fraction  fields  of  the  ethylene 
flame  simulations  with  and  without  the  EDC  model  are  shown  in  Figure  46.  The  effect  of 
turbulence-chemistry  interaction  on  the  ignition  of  the  flame  is  even  more  significant  for 
the  ethylene  flame.  For  the  case  with  laminar  chemistry,  the  flame  was  ignited  half  way 
through  the  computational  domain,  while  for  EDC  model,  the  flame  is  barely  started  at 
the  end  of  the  computational  domain.  The  distribution  of  C2H4  in  both  cases  is  similar 
because  the  amount  of  C2H4  consumed  in  the  two  cases  is  very  small  due  to  the  high 
speed  of  the  jet. 
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Figure  46.  Temperature  and  C2H4  Mass  Fraction  Fields  Predicted  with  EDC  Model 
and  Laminar  Chemistry. 


The  CPU  timings  of  the  EDC  model  for  the  hydrogen  and  ethylene  flames  are  shown  in 
Table  9. 
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Table  9.  CPU  Time  (in  seconds)  Per  Time  Step  for  Simulations  with  EDC  Model 
Compared  with  Laminar  Chemistry. 


Mechanisms 

Lam.  Chem. 

w/o  ISAT 

Lam.  Chem. 

w/  ISAT 

EDC 

w/o  ISAT 

EDC 

w/  ISAT 

H2,  7-specie 

0.254 

0.147 

11.2 

0.431 

C2H4,  20-specie 

18.7 

0.51 

908.1 

3.21 

As  previously  discussed,  there  are  only  two  steady-state  species  in  the  reduced  hydrogen 
mechanism.  Calculation  with  and  without  ISAT  is  relatively  inexpensive.  The  CPU  time 
including  the  EDC  model  without  ISAT  is  1 1.2  seconds  per  time  step,  44  times  that  for 
laminar  chemistry  without  ISAT.  With  ISAT,  the  CPU  time  for  the  EDC  model  is 
reduced  to  0.43 1  seconds,  a  speedup  factor  of  26. 

The  20-specie  C2H4  mechanism  is  significantly  more  expensive.  The  CPU  time  of  the 
EDC  model  without  ISAT  is  908.1  seconds  per  iteration.  It  is  reduced  to  3.21  seconds  per 
time  step  with  ISAT.  The  speedup  factor  is  282.  The  CPU  cost  for  a  simulation  with  the 
EDC  model  and  ISAT  for  the  ethylene  flame  is  approximately  6  times  that  using  laminar 
chemistry  with  ISAT.  The  tremendous  CPU  savings  from  the  ISAT  makes  it  feasible  to 
perform  a  3-D  simulation  with  a  practical  fuel. 

9.2.2  PDF  Post-Processor 

The  axisymmetric  hydrogen  flame  was  chosen  as  the  test  case  for  the  PDF  post-processor 
to  limit  the  tremendous  CPU  cost  involved  in  PDF  simulations.  The  reduced  chemical 
mechanism  is  the  seven-specie  H2  mechanism.  A  curvilinear  coordinate  system  is  used  in 
VULCAN.  However,  the  PDF  post-processor  requires  Cartesian  coordinates.  The 
predicted  velocity,  density,  and  turbulent  kinetic  energy  and  its  dissipation  rate  from  the 
VULCAN  simulation  were  interpolated  to  a  61  by  51  Cartesian  grid.  Those  variables 
were  fixed  in  the  post-processor,  while  the  species  concentrations  and  temperature  were 
recalculated.  Two  mixing  models  were  applied  in  the  post-processor:  a  well-mixed  model 
(no  turbulence-chemistry  interaction)  and  a  modified  Curl’s  model  (Curl,  1963).  For 
Curl’s  model,  the  model  constant  C  ,  which  is  the  mechanical-to-scalar  timescale  ratio,  is 

set  to  2,  a  typical  value  in  PDF  simulations. 

The  predicted  mean  temperature  and  H2  and  02  mass  fractions  with  the  well-mixed 
model  are  shown  in  Figure  47.  The  results  agree  well  with  the  results  from  VULCAN 
using  laminar  chemistry  (see  Figure  41),  although  the  PDF  post-processor  predicts  a 
slightly  larger  flame  liftoff  height.  Since  different  numerical  schemes  are  employed  in 
VULCAN  and  the  post-processor,  a  perfect  match  is  not  expected. 
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Figure  47.  Contour  Plots  of  Mean  Temperature  and  H2  and  02  Mass  Fractions  for 
the  Well-Mixed  Model. 


The  results  predicted  with  the  modified  Curl’s  model  are  shown  in  Figure  48.  It  is 
surprising  to  see  that  the  flame  is  quenched  in  this  case.  The  maximum  temperature  in  the 
jet  was  slightly  higher  than  that  of  the  air  stream.  To  understand  the  reason  for  the 
extinction,  the  turbulent  time  scale,  t ,  which  is  calculated  from  t  k  /  ,  is  plotted  in 

Figure  49.  This  variable  was  used  to  determine  the  mixing  frequency  in  the  Curl’s  model, 
which  is  given  by 


H2  mass  fraction  *0.0 


02  mass  fraction 


0.24 


0-00 


Figure  48.  Contour  Plots  of  Mean  Temperature  and  H2  and  02  Mass  Fractions  for 
the  Modified  Curl’s  Model. 


Figure  49.  The  Turbulent  Time  Scale,  t  k  /  (in  seconds). 
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Inside  the  jet  shear  layer  where  mixing  between  the  hydrogen  and  air  takes  place,  the 
time  scale  is  on  the  order  of  10'2  seconds.  The  residence  time,  defined  as  the  time  for  a 
fluid  particle  to  walk  through  the  whole  computational  domain  from  the  inlet  to  the 
outlet,  is  approximately  1.5  xlO"4  seconds  (the  length  of  the  computational  domain  is 
0.28  m,  and  the  averaged  axial  velocity  is  approximately  2000  m/s).  Thus,  residence  time 
is  2  orders  of  magnitude  smaller  than  t .  This  implies  that  the  molecular-level  mixing 

between  the  particles  from  the  hydrogen  and  air  streams  is  very  slow  compared  to  the 
convection  speed  of  the  jet. 


We  define  a  new  variable,  unmixedness: 

1  //’ 


(15) 


where  /  and  /  are  the  mean  and  variation  of  the  mixture  fraction,  respectively.  The 
variable  is  0  in  the  well-mixed  region,  and  it  is  1  in  the  region  where  the  fluids  are  totally 
unmixed.  Figure  50  shows  the  unmixedness  for  the  well-mixed  model  and  modified 
Curl’s  model.  For  the  well-mixed  model,  is  0  everywhere.  For  the  modified  Curl’s 
model,  it  is  1  inside  the  jet  shear  layer,  and  0  in  the  pure  air  and  hydrogen  regions. 


Figure  50.  Contour  Plots  of  Unmixedness  for  Well-Mixed  Model  (top)  and  Modified 
Curl’s  Model  (bottom). 

The  turbulent  time  scale  defined  by  kl  as  predicted  in  VULCAN,  is  large.  It  results  in 
the  extinction  of  the  hydrogen  flame.  To  study  the  effect  of  the  time  scale  on  the  PDF 
calculation,  the  time  scale  was  artificially  decreased  by  a  factor  of  100  such  that  its  value 
inside  the  jet  shear  layer  is  in  the  same  order  of  magnitude  as  the  residence  time.  The 
predicted  flame  using  the  modified  Curl’s  model  is  shown  in  Figure  51.  With  the 
improved  molecular-level  mixing  between  the  particles,  the  flame  is  stabilized.  The 
liftoff  height  is  larger  than  that  of  the  well-mixed  case,  which  is  consistent  with  the 
predictions  of  the  EDC  model. 
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Figure  51.  Predicted  Mean  Temperature  and  H2  and  02  Mass  Fractions  Using  an 
Artificially  Increased  Turbulent  Time  Scale. 


These  simulations  indicate  that  the  turbulent  time  scale  is  very  important  in  determining 
the  molecular-level  mixing  process  in  the  PDF  methods. 

9.2.3  Difference  Between  Turbulence-Chemistry  Interaction  Models 

Evans  et  al.  (1978,  Appendix  B)  described  a  method  for  calculating  the  profdes  of  k  and 
based  on  Prandtl’s  mixing  length  equation.  Although  the  method  was  applied  only  to 
provide  an  initial  profde  for  their  simulations,  it  can  also  give  us  an  estimation  of  the 
magnitudes  of  k  and  throughout  a  jet  shear  layer.  The  turbulent  kinetic  energy  and  its 
dissipation  rate  estimated  by  Evans  et  al.  peaked  at  approximately  105  m2/s2  and  1011 
m2/s3,  respectively  (See  Evans  et  al.,  1978,  Figure  1  (b)).  While  the  peak  value  of 
k  predicted  by  VULCAN  is  in  the  same  order  of  magnitude  as  the  Evans  estimate,  the 
value  of  is  significantly  lower,  with  peak  values  of  only  107  m2/s3.  The  apparent  under¬ 
prediction  of  in  VULCAN  results  in  a  very  low  value  of  mixing  frequency,  which 
results  in  the  flame  being  quenched  in  the  PDF  post-processor  simulations. 

In  the  EDC  model,  k  and  are  used  to  calculate  the  volume  fraction  of  flow  occupied 
by  fine  structures,  ,  and  the  lifetime  of  the  fine  structures,  .  For  low  values  of  ,  and 
approach  1  and  0,  respectively.  From  Eq.  (7-6),  the  turbulent  reaction  rate  approaches 
the  laminar  reaction  rate,  which  is  consistent  with  the  results  shown  earlier  in  this  section. 
It  is  expected  that,  with  a  higher  calculated  value  of  ,  the  flame  structures  predicted 
using  EDC  model  and  PDF  method  will  be  much  closer  to  each  another. 

In  summary,  simulations  with  both  the  EDC  and  PDF  methods  showed  the  effect  of 
turbulent  mixing  on  flame  structure,  i.e.  larger  lift-off  height  and  a  lower  maximum 
temperature  inside  the  flame.  Both  models  rely  on  the  accuracy  of  the  turbulence  model 
within  the  CFD  code  to  predict  the  turbulent  kinetic  energy  and  dissipation  rate. 
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10.0  3-D  Scramjet  Simulations 

Two  of  the  previously  described  reduced  mechanisms,  those  for  ethylene  and  ^-heptane 
have  been  used  to  run  3-D  simulations  of  a  scramjet  cavity  flameholder.  Results  of  these 
simulations  are  reported  in  this  section.  These  simulations  were  performed  on  a  Compaq 
ES-45  computer  using  sixteen  processors.  Each  simulation  required  about  250,000 
iterations  and  about  one  month  of  computer  time  to  converge.  The  grid  for  these 
simulations  had  334,558  cells. 

10.1  Ethylene 

The  ethylene  reduced  chemical  kinetic  mechanism  previously  incorporated  into 
VULCAN  was  tested  for  a  3-D  simulation.  The  test  model  is  a  cavity-based  flameholder 
with  an  injection  port  at  the  back  of  the  cavity-ramp  wall,  which  corresponds  to  the 
experimental  research  of  AFRL/PRA  Test  Cell  19  facility.  Although  there  were  10 
injection  ports  in  the  back  cavity-ramp  wall,  current  numerical  simulation  consisted  of 
the  study  of  the  flow  fields  of  one  of  the  injector  ports,  as  shown  in  Figure  52.  In 
addition,  no  side-wall  effects  are  considered  in  these  simulations.  The  freestream  Mach 
number  is  2.0,  with  the  total  temperature  and  pressure  of  ethylene  fuel  at  the  injector  exit 
at  31 1  K  and  12  psia,  respectively. 

In  the  numerical  simulation,  the  ethylene  fuel  was  first  injected  into  the  cavity  so  that  the 
mixing  between  the  fuel  and  air  in  the  cavity  could  be  established.  An  ignition  source  was 
then  applied  in  the  cavity  region  to  initiate  combustion.  Once  combustion  had  been 
established  based  on  increased  temperature  and  production  of  CO2  and  H2O  in  the  cavity 
region,  the  ignition  source  was  turned  off.  The  simulation  was  then  continued  to 
determine  whether  the  flame  in  the  cavity  would  be  self-sustaining.  In  the  Test  Cell  19 
experimental  setup,  spark  plugs  were  used  for  ignition.  Experimental  results  have  shown 
that  the  flame  is  able  to  self-sustain  under  these  conditions.  Numerical  results  using  the 
CARM  reduced  mechanism  also  showed  that  the  flame  was  able  to  sustain  itself  in  the 
cavity,  as  shown  in  Figures  53  through  57.  Previous  efforts  to  simulate  these  conditions 
with  global  ethylene  kinetics  models  have  incorrectly  predicted  flame  extinction  at  these 
conditions.  This  demonstrates  CARM’s  ability  to  incorporate  more  chemical  kinetic 
information  into  a  model  than  curve  fitting  of  global  reaction  rates  can. 

The  temperature  contours  (Figure  53)  show  high  temperature  regions  in  the  cavity,  which 
indicate  combustion  taking  place.  Note  that  the  low  temperature  zone  in  the  center  of  the 
cavity  is  due  to  very  fuel  rich  conditions.  This  can  be  observed  from  the  E.R.  contours 
(Figure  54).  CO2  and  OH  mass  fractions  are  depicted  in  Figures  55  and  56.  Regions  of 
high  concentrations  of  these  species  show  that  the  chemical  reaction  is  taking  place  near 
the  stoichiometric  surface  (shown  in  Figure  54).  The  pressure  contours  and  streamlines 
are  shown  in  Figure  57.  The  recirculation  region  in  the  cavity  facilitates  mixing  of  fuel 
and  air  as  well  as  hot,  radical-containing  flame  gases  with  unbumed  fuel-air  mixture, 
which  sustains  the  flame,  even  under  high  Mach  number  freestream  conditions. 
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Figure  52.  Computational  Domain  with  Cavity  and  Injector  Port. 
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Figure  53.  Temperature  Contours  (K)  in  the  Vicinity  of  the  Cavity. 
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Figure  54.  E.R.  Contours  in  the  Cavity  Region. 
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Figure  55.  CO2  Mass  Fraction  Contours  in  the  Cavity  Region. 
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Figure  56.  OH  Mass  Fraction  Contours  in  the  Cavity  Region. 
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Figure  57.  Pressure  Contours  and  Streamlines  in  the  Cavity  Region. 

10.2  7V-Heptane 

The  /7-heptane  reduced  mechanism  described  in  previous  reports  has  been  used  to  model 
a  3-D  flameholder  case.  The  conditions  were  similar  to  those  of  the  ethylene  case.  Once 
again  stability  problems  due  to  the  reduced  mechanism  were  not  observed,  and  a  self- 
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sustaining  simulated  flame  was  achieved.  Results  of  this  calculation  are  shown  in 
Figures  58  through  61. 

The  test  model  is  a  cavity-based  flameholder  with  an  injection  port  at  the  back  of  the 
cavity-ramp  wall,  similar  to  the  ethylene  case.  However,  no  experimental  data  is 
available  for  comparison  for  the  /7-heptane  case.  In  the  ethylene  simulation,  the  total 
pressure  and  temperature  of  the  injector  was  provided  by  AFRL/PRA  Test  Cell  19 
facility,  and  the  data  corresponds  to  82.74  kPa  (12  psia)  and  316.46  K  (110  °F), 
respectively. 

Extensive  effort  was  required  to  choose  appropriate  boundary  conditions  at  injector  port 
entrance  for  the  «-heptane  gas.  Two  analytical  approaches  were  used  to  calculate  the 
pressure  coming  out  of  the  injector,  while  keeping  the  total  temperature  (316.46  K)  the 
same  as  in  the  ethylene  case.  The  first  attempt  was  to  assume  that  the  mass  flow  rate  and 
injector  area  of  the  ethylene  and  heptane  gases  are  the  same  out  of  the  injector  port.  Then 
the  total  pressure  of  the  //-heptane  gas  can  be  computed  using  Equation  (16),  where  w  and 
A  are  the  mass  flow  rate  and  injector  area,  respectively: 

The  values  of  and  R  were  1.26  and  83.1493  J/kg  K,  based  on  the  temperature  and  the 
molecular  weight  of  the  //-heptane  gas.  The  total  pressure  was  found  to  be  43.33  kPa. 
These  total  properties  were  used  at  the  entrance  of  the  injector  boundary  condition. 
However,  in  this  simulation,  the  static  pressure  at  the  injector  entrance  was  lower  than  at 
the  injector  exit,  i.e.,  at  the  flameholding  cavity.  Thus,  the  air  inside  the  cavity  was 
flowing  into  the  injector  port,  instead  of  the  heptane  gas  flowing  into  the  cavity. 


The  second  approach  used  the  isentropic  relation  and  assumed  that  the  flow  was  choked 
at  the  injector  exit.  The  static  pressure  was  set  to  be  56  kPa,  which  corresponded  to  the 
injector  exit  for  the  ethylene  case.  Thus,  the  total  pressure  was  found  to  be  101.3  kPa, 
while  the  total  temperature  remained  at  316.46  K.  In  this  case,  the  numerical  solution 
showed  that  the  flow  was  choked  inside  the  injector  port.  The  //-heptane  gas  was  initially 
at  a  supersonic  flow  and  became  subsonic  inside  injector  port.  Eventually,  this  solution 
diverged.  Different  total  pressure  values  were  tested  for  the  injector  entrance. 
Unfortunately,  these  solutions  were  either  not  promising  and/or  diverging. 

To  reduce  the  complexity  of  the  fuel  inflow  boundary  condition,  the  injector  itself  was 
eliminated  from  the  model,  and  a  fixed  boundary  condition  was  imposed  at  the  injector 
exit.  A  top-hat  velocity  profile  was  assumed.  Initially,  the  total  pressure  was  set  to  43.33 
kPa,  which  provided  the  same  mass  flow  rate  as  the  ethylene  case.  However,  the  solution 
for  these  conditions  did  not  converge  due  to  the  low  static  pressure  at  the  vicinity  of  the 
injector  exit.  After  trial-and-error,  the  static  pressure  was  increased  to  160  kPa  with  exit 
velocity  of  1 80  m/s.  The  simulation  using  this  boundary  condition  has  shown  promise. 
First,  a  nonreacting  flow  simulation  was  performed  to  establish  a  mixing  solution.  Then, 
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the  ignition  zone  was  increased  gradually  inside  the  cavity  with  a  constant  static 
temperature  of  1500  K  to  start  a  reacting  solution.  Once  combustion  had  been  established 
based  on  increased  temperature  and  production  of  CO2  in  the  cavity  region,  the  ignition 
source  was  turned  off  to  determine  whether  the  flame  in  the  cavity  would  self-sustain. 
Although  the  computation  has  not  completed,  the  preliminary  result  shows  that  the  flame 
was  able  to  self-sustain  in  the  cavity,  as  shown  in  Figures  58  through  61. 

The  temperature  contours  (Figure  58)  show  high  temperature  regions  in  the  cavity,  which 
indicate  that  combustion  is  taking  place  there.  The  high  temperature  zones  are  located 
near  the  front  of  the  cavity  and  center  of  the  bottom  walls.  The  temperature  is  lower  at  the 
center  of  the  cavity  due  to  fuel-rich  conditions.  This  can  be  observed  from  the  E.R. 
contours  (Figure  59)  Mass  fractions  of  CO2  and  H  are  depicted  in  Figures  60  and  61. 
Regions  with  high  concentrations  of  these  species  show  where  chemical  reactions  were 
taking  place.  This  simulation  demonstrates  the  capability  of  CARM  to  produce  a  reduced 
mechanism  for  a  large  hydrocarbon  fuel  that  can  be  implemented  and  run  in  a  3-D 
VULCAN  case. 
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Figure  58.  Temperature  Contours  (K)  in  the  Vicinity  of  the  Cavity. 
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Figure  59.  E.R.  Contours  in  the  Cavity  Region. 
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Figure  60.  CO2  Mass  Fraction  Contours  in  the  Cavity  Region. 
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Figure  61.  H  Mass  Fraction  Contours  in  the  Cavity  Region. 
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11.0  Conclusions 

This  project  has  addressed  existing  limitations  regarding  our  ability  to  represent  complex 
hydrocarbon  combustion  chemistry  within  3-D  simulations  of  high-speed  combustion  of 
complex  aviation  fuels.  This  project  focused  on  the  implementation  of  several  key  recent 
developments  into  existing  state-of-the-art  CFD  software  that  is  currently  utilized  in  the 
design  of  scramjet  combustors.  The  developments  that  were  utilized  in  this  project  are  as 
follows: 

Detailed  chemistry  for  aviation  fuels 

Software  (CARM)  for  generating  accurate  reduced  mechanisms  for  aviation  fuels 
Sophisticated  tabulation  techniques  (ISAT)  that  reduce  the  CPU  time  required  for 
modeling  hydrocarbon  combustion 

Techniques  (EDC  &  PDF)  to  model  the  impacts  of  turbulence  on  calculation  of 
the  mean  reaction  rates 

Software  (VULCAN)  for  modeling  high-speed  hydrocarbon  combustion. 

This  project  integrated  these  components  to  demonstrate  a  process  for:  1)  developing 
accurate  reduced  chemical  kinetic  mechanisms  for  aviation  fuels,  2)  integrating  them 
with  ISAT  and  EDC  into  VULCAN,  and  3)  simulation  of  combusting  flows  with 
complex  aviation  fuel  chemistry  in  scramjet  combustor  geometries.  Data  for  the  2-D  and 
3-D  cases  is  very  limited  particularly  for  aviation  fuels.  This  project  focused  on  utilizing 
simple  configurations  for  verification  of  accuracy  in  cases  where  data  were  available. 

The  key  results  from  the  Phase  II  work  included: 

Development  of  accurate  reduced  mechanisms  for  combustion  kinetics  of 
hydrogen,  ethylene,  n-heptane,  and  a  six-specie  surrogate  blend  representing  JP-8 
Improvements  to  CARM  software  for  automated  reduction  of  chemical  kinetic 
mechanisms  resulting  in  improved  robustness  and  decreased  computational  times 
Implementation  of  four  reduced  mechanisms  into  VULCAN 
Implementation  of  ISAT  chemical  rate  tabulation  technique  into  VULCAN 
Implementation  of  EDC  turbulence  model  (with  ISAT)  into  VULCAN 
Completion  of  simulations  demonstrating  improved  chemical  descriptions  from 
reduced  mechanisms  compared  to  previously  used  global  chemistry  models 
ISAT  results  showing  speed  increases  up  to  a  factor  of  42 

2- D  combustion  simulations  in  VULCAN  demonstrating  this  approach  for 
hydrogen,  ethylene,  n -heptane,  and  JP-8  (up  to  35  species  in  the  reduced 
mechanism) 

3- D  scramjet  combustor  simulations  run  for  ethylene  and  « -heptane  (up  to  20 
species  in  the  reduced  mechanism),  demonstrating  that  this  process  can  be  used 
for  CFD  simulations  of  realistic  3-D  scramjet  geometries  with  aviation  fuels. 

The  CARM  software  has  proven  to  be  a  powerful  tool  for  reducing  large  detailed 
kinetic  mechanisms  modeling  the  combustion  of  hydrocarbon  fuels  from  the  literature 
to  a  form  that  can  be  implemented  into  a  CFD  code.  ISAT  is  a  necessary  companion 
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technology  which  is  needed  to  reduce  the  computational  cost  of  source  term 
evaluation  when  using  CARM-reduced  mechanisms. 

Modeling  turbulence-chemistry  interactions  for  high-speed  reacting  flows  remains  a 
challenge  due  to  both  the  computational  cost  of  the  models  and  the  uncertainties  of 
the  models  themselves.  Both  the  EDC  and  PDF  models  tested  during  this  project 
gave  the  expected  results  of  increasing  the  flame  liftoff  distance  and  reducing  the 
maximum  flame  temperature.  However,  the  models  gave  very  different  predictions. 
These  differences  can  be  traced  to  the  ways  the  models  use  the  information,  such  as 
turbulent  time  scales,  from  the  fluid  turbulence  model.  The  results  suggest  that 
further  work  is  needed  to  improve  both  turbulence  models  and  turbulence-chemistry 
interaction  models  for  high-speed  reacting  flows. 
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12.0  Future  Directions 

REI  has  been  actively  developing  and  applying  reduced  chemistry  to  3-D  simulation  for  a 
wide  range  of  industrial  applications  involving  low-speed  combustion  and  incineration 
and  will  continue  to  do  so.  There  are  also  a  number  of  potential  follow-on  opportunities 
associated  with  simulation  of  high-speed  reacting  flows  related  to  the  work  completed  on 
this  project  that  could  be  pursued.  Some  specific  areas  are  summarized  here. 

12.1  Development  of  a  Reduced  Mechanism  for  JP-7 

JP-7  is  a  highly  refined  hydrocarbon  fuel  originally  developed  for  the  SR-7 1  (Edwards  & 
Maurice  2001)  but  now  of  interest  as  a  scramjet  fuel.  Kinetics  for  JP-7  have  been 
developed  (Mawid  et  al.,  2003)  but  they  are  proprietary  and  their  distribution  may  be 
problematic.  It  is  likely  that  a  suitable  surrogate  could  be  formulated  using  compounds 
available  in  the  Violi  et  al.  (2002)  mechanism,  from  which  a  JP-8  reduced  mechanism 
was  successfully  created  during  this  project.  It  is  also  possible  that  the  present  JP-8 
reduced  mechanism  could  work  for  JP-7  with  a  surrogate  fuel  containing  species  in  the 
JP-8  mechanism  but  in  different  proportions.  Agreement  among  interested  parties  on  a 
suitable  detailed  mechanism  and  fuel  surrogate  would  be  necessary  before  starting  work 
on  the  development  of  a  reduced  mechanism. 

12.2  Incorporation  of  ISAT/CARM  Methodology  into  Metacomp’s  CFD++ 

VULCAN  is  not  the  only  CFD  code  used  by  Air  Force  personnel  for  simulation  of  high¬ 
speed  combusting  flows.  Metacomp’s  CFD++  is  a  well-developed  commercial  code  that 
could  see  significant  benefit  from  incorporation  of  IS  AT  and  reduced  mechanism 
capability.  This  would  require  an  agreement  between  REI  and  Metacomp,  possibly 
allowing  REI  access  to  at  least  some  source  code.  It  may  also  require  a  C  code  output 
option  from  CARM. 

12.3  Improvements  in  ISAT  Table  Access  Methods 

Several  possibilities  exist  for  improving  the  efficiency  of  the  ISAT  algorithm.  Among 
these  are  periodic  table  restructuring  and  multitable  methods.  Possibilities  also  exist  for 
improving  the  performance  of  ISAT  on  multiple  processors.  Investigation  into  the 
sensitivity  of  accuracy  and  speedup  to  alternate  definitions  of  the  Region  of  Accuracy 
could  also  be  beneficial. 

12.4  Turbulence-Chemistry  Interaction  Methods  for  High-Speed  Flows 

This  is  a  difficult  but  important  area  of  research.  New  ideas  that  improve  on  a  constant 
Schmidt  number  approximation  would  be  valuable  to  the  high-speed  flow  community. 
Novel  ideas  could  involve  aspects  of  existing  models  for  turbulent  mixing,  including 
Kerstein’s  LEM. 
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List  of  Acronyms 


Acronym 

Description 

2-D 

two  dimensional 

3-D 

three  dimensional 

BSL 

baseline 

CARM 

computer-assisted  reduction  mechanism 

CFD 

computational  fluid  dynamic  (or  dynamics) 

CPU 

central  processing  unit 

DDASAC 

Double  precision  Differential-Algebraic  Sensitivity  Analysis  Code 

DoD 

Department  of  Defense 

DR 

detailed  reduction 

E.R. 

equivalence  ratio 

EBU 

eddy  breakup 

EDC 

eddy  dissipation  concept 

EOA 

ellipsoid  of  accuracy 

GRI 

Gas  Research  Institute 

H-0 

hydrogen-oxygen 

ISAT 

in  situ  adaptive  tabulation 

KIVA 

not  an  acronym 

LEM 

linear  eddy  model 

LU 

lower-upper 

NSF 

National  Science  Foundation 

ODE 

ordinary  differential  equation 

PDF 

probability  density  function 

PFR 

plug  flow  reactors 

PSE 

problem  solving  environment 

PSR 

perfectly  stirred  reactor 

r/D 

radius  divided  by  diameter 

QSS 

quasi-steady  state 

REI 

Reaction  Engineering  International 

UNICORN 

UNsteady  Ignition  and  COmbustion  with  ReactioNs 

VULCAN 

Viscous  Upwind  aLgorithm  for  Complex  flow  ANalysis 
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